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We introduce numerical methods for simulating the diffusive motion of rigid bodies of 
arbitrary shape immersed in a viscous fluid. We parameterize the orientation of the bodies 
using normalized quaternions, which are numerically robust, space efficient, and easy to 
accumulate. We construct a system of overdamped Langevin equations in the quaternion 
representation that accounts for hydrodynamic effects, preserves the unit-norm constraint on 
the quaternion, and is time reversible with respect to the Gibbs-Boltzmann distribution at 
equilibrium. We introduce two schemes for temporal integration of the overdamped Langevin 
equations of motion, one based on the Fixrnan midpoint method and the other based on a 
random finite difference approach, both of which ensure the correct stochastic drift term is 
captured in a computationally efficient way. We study several examples of rigid colloidal 
particles diffusing near a no-slip boundary, and demonstrate the importance of the choice of 
tracking point on the measured translational mean square displacement (MSD). We examine 
the average short-time as well as the long-time quasi-two-dimensional diffusion coefficient of 
a rigid particle sedimented near a bottom wall due to gravity. For several particle shapes 
we find a choice of tracking point that makes the MSD essentially linear with time, allowing 
us to estimate the long-time diffusion coefficient efficiently using a Monte Carlo method. 
However, in general such a special choice of tracking point does not exist, and numerical 
techniques for simulating long trajectories, such as the ones we introduce here, are necessary 
to study diffusion on long timescales. 


I. INTRODUCTION 

The Brownian motion of rigid bodies suspended in a viscous solvent is one of the oldest subjects 
in nonequilibrium statistical mechanics, and is of crucial importance in a number of applications 
in chemical engineering and materials science. Examples include the dynamics of passive m 
or active pMo] particles in suspension, the dynamics of biomolecules in solution mHEg, the 
design of novel nano-colloidal materials m, and others. At the mesoscopic scales of interest, the 
erratic motion of individual molecules in the solvent drives the diffusive motion of the suspended 
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particles. The number of degrees of freedom necessary to simulate this motion directly using 
Molecular Dynamics (MD) is large enough to make this approach prohibitively expensive. Instead, 
the Brownian dynamics approach captures the effect of the solvent through a mobility operator, 
and thermal fluctuations are modeled using appropriate stochastic forcing terms. In previous 
work P!> we used a computational fluid solver and immersed boundary techniques to simulate 
the diffusive motion of spherical particles including hydrodynamic interactions. The fluctuating 
immersed boundary method developed in US! is suitable for minimally-resolved computations in 
which only the translational degrees of freedom are kept and hydrodynamics is resolved at a far- 
field level assuming the particles are spherical. Novel methods are, however, required to model 
the behavior of particles with nontrivial shapes such as rigidly-fused colloidal clusters El SI or 
colloidal boomerangs [I]. In this paper, we show how to include rotational degrees of freedom 
in the overdamped Langevin equations of motion for rigid bodies suspended in a viscous fluid, 
develop specialized temporal integrators for these equations, and apply them to a number of model 
problems. 

One of the important goals of our work is to develop an overdamped formulation and associated 
numerical algorithms that apply when the hydrodynamic mobility (equivalently, resistance) de¬ 
pends strongly on the configuration. Many previous works have focused on the rotational diffusion 
of a single isolated rigid body in an unbounded domain. However, in practice, rigid particles diffuse 
either in a suspension, in which case they interact hydrodynamically with other particles, or near 
a boundary such as a microscope slide or the walls of a slit channel, in which case they interact 
hydrodynamically with the boundaries. Here we consider a general case of a rigid body performing 
translational and rotational Brownian motion in a confined system, specifically, we numerically 
study particles sedimented close to a single no-slip boundary. This is of particular relevance to 
recent experimental studies of the diffusive motion of colloidal particles that are much denser than 
water and thus sediment close to the microscope slide (glass plate) [THS]. 

When writing the equations of motion for a rigid body one must first choose how to represent 
the orientation of the body. For bodies with a high degree of symmetry one can use simple 
representations of orientation, for example, for axisynunetric particles (e.g., rigid rods) in three 
dimensions one can use two polar angles or a unit vector to represent the orientation of the axis of 
symmetry [5, 16,(TJ|. More complex (biaxial or skewed) particle shapes |H[6], or asymmetrically 
patterned particles of symmetric shapes |3], as common in active particle suspensions mi , require 
describing the complete orientation of the rigid bodies. Mathematically, the orientation of a general 
rigid body in three dimensions is an element of the rotation group 50(3); the group of unitary 3x3 
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matrices of unit determinant (rotation matrices). This group can be parameterized in a number of 
ways, the most fundamental one representing elements of this group by an orientated rotation angle, 
represented as a three-dimensional vector 4 >, the direction of which gives an axes of rotation relative 
to a reference configuration, and the magnitude of which gives an angle of rotation around that axes. 
Prior work on rotational Brownian motion in the overdamped regime has considered the use of Euler 
angles mm, oriented rotation angles m, as well as a number of other representations USED]. 
Each of these representations has its own set of problems, notably, most of them have singularities 
or redundancies (which can be avoided in principle with sufficient care), lead to complex analytical 
expressions involving potentially expensive-to-evaluate trigonometric functions, or require a large 
amount of storage (e.g., a rotation matrix with 9 elements). Furthermore, with the exception of 
P33 HI 119) . most prior work on rotational diffusion either assumes that the mobility does not 
depend on configuration [21), focuses on cases where tracking a single axes is sufficient to describe 
the Brownian motion mmm, or is not careful in handling the stochastic drift terms necessary 
when the rotational mobility is dependent on the position and orientation of the body. 

In molecular dynamics circles [231425] . it is well-known that a robust and efficient representation 
of orientation is provided by unit quaternions, which are unit vectors in four dimensions (i.e., points 
on the unit 4-sphere). This representation contains one redundant degree of freedom (four instead of 
the minimal required of three), however, it is free of singularities and thus numerically robust, and, 
as we will see, leads to a straightforward formulation that is simple to work with both analytically 
and numerically. In some sense, the quaternion representation is a direct generalization to bi-axial 
bodies of the standard representation used in Brownian Dynamics of uni-axial particles [5I|, namely, 
a unit vector in three dimensions. That common representation is also redundant (only two polar 
angles are required to describe a direction in three dimensions), however, it offers many advantages 
over more compressed representations such as polar angles, and is thus the representation of choice. 
Following the submission of this manuscript, we learned of a very recent work by Hie et al that also 
uses quaternions in an overdamped Langevin equation for the motion of a general rigid body in 
bulk [26]; earlier work m has also used quaternions but without carefully considering the required 
stochastic drift terms. 

We consider the overdamped regime, where the timescale of momentum diffusion in the fluid is 
much shorter than the timescale of the motion of the rigid bodies themselves. Formally, this regime 
corresponds to the limit of infinite Schmidt number [23]. Neglecting inertia, we track only the 
positions and orientations of the immersed bodies, deriving evolution equations for the quaternion 
representation. This Langevin system exhibits the correct deterministic dynamics and preserves the 
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Gibbs-Boltzmann distribution in equilibrium, properly restricted to the unit quaternion 4-sphere. 
Integrating these equations proves challenging primarily due to the presence of the stochastic drift 
term that arises from the configuration-dependent mobility; this issue is identified theoretically in 
Appendix C in [26] but that work is focused on unconfined particles for which a key stochastic 
drift term vanishes (see (C21) in [25]). The standard approach to handling the stochastic drift 
term is Fixman’s method, requiring a costly application of the inverse of the mobility which in 
some cases is not directly computable. As an alternative, we employ a recently-proposed Random 
Finite Difference (RFD) scheme |151129] for approximating the drift; this approach only requires 
application of the mobility and its “square root” but not the inverse of the mobility. 

We perform a number of numerical experiments in which we simulate the Brownian motion of 
rigid particles sedimented near a wall in the presence of gravity, as inspired by recent experimental 
studies of the diffusion of asymmetric spheres [3], clusters of spheres im and boomerang colloids 
mi- In the first example, we study a tetramer formed by rigidly connecting four colloidal spheres 
placed at the vertices of a tetrahedron, modeling colloidal clusters that have been manufactured in 
the lab In the second example, we study the rotational and translational diffusion of 

an asymmetric colloidal sphere with center of mass displaced from the geometric center, modeling 
recently-manufactured “colloidal surfers” [8] in which a dense hematite cube is embedded in a 
polymeric spherical particle. In the last example we study the quasi two-dimensional diffusive 
motion of a dense boomerang colloid sedimented near a no-slip boundary, as inspired by recent 
experiments P 0 [32]. We computationally demonstrate the crucial importance of the choice 
of tracking point when computing the translational diffusion coefficient. In particular, we show 
that with a suitable choice of the origin around which torques are expressed, one can obtain 
an approximate but relatively accurate formula for the effective long-time diffusion coefficient in 
the directions parallel to the boundary. However, we are unable to reach a precise and definite 
conclusion about the optimal choice of tracking point even for quasi-two-dimensional diffusion, 
since for all shapes studied here and in existing experiments the center of hydrodynamic stress 
and the center of mobility are too close to each other to be distinguished. In the more general 
case, our results indicate that there is no exact closed-form expression for the long-time quasi-two- 
dimensional coefficient, and numerical methods for simulating trajectories are necessary in order to 
study the long-time diffusive dynamics of even a single rigid body in the presence of confinement. 

This paper is organized as follows. In Section [TT] we formulate the equations of motion for 
rigid bodies with translation and rotation, giving a brief background on the use of quaternions to 
parameterize orientation. Section |ITT| introduces temporal integrators for these equations, including 
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a Fixman scheme, as well as a RFD scheme that approximates the stochastic drift using only 


applications of the mobility. We perform numerical tests of our schemes in Section IV to verify 
that we can correctly simulate the dynamics of a rigid body near a no-slip boundary, and study 
the influence of the choice of tracking point on the MSD. Finally, we give concluding thoughts and 
discuss future directions in Section IVl Technical details are handled in Appendices. 


II. LANGEVIN EQUATIONS FOR RIGID BODIES 

In this section, we formulate Langevin equations for rigid bodies performing rotational and 
translational diffusion. We begin by formulating an overdamped Langevin equation for rotational 
diffusion using a unit quaternion representation of rigid-body orientation. For the remainder of this 
section, we will assume that we know how to compute the configuration dependent hydrodynamic 
mobilities needed for our equations. These mobility matrices are applied to vectors of forces and 
torques to compute the resulting linear and angular velocities of the immersed rigid bodies. In 
future work, we will develop algorithms for computing these objects on the fly using a computational 
fluid solver as in the Fluctuating Immersed Boundary (FIB) method [13], as we discuss in more 
detail in Section [Vj 

Our goal is to formulate an equation for the evolution of the orientation of a rigid body. It is 
important that the resulting system has the correct deterministic term, that it is time reversible 
with respect to the correct Gibbs-Boltzmann distribution in equilibrium, and that it preserves the 
constraint that the quaternion has unit norm. Before we accomplish this goal, we briefly review 
some required facts about quaternions. 


A. Quaternions 

Describing the orientation of a rigid body in three dimensions can be done in many ways. 
Rotation matrices are perhaps the most straightforward approach to accomplish this task, but they 
require the use of 9 floating point numbers to parameterize a 3 dimensional space. Additionally, 
accumulation of numerical errors over many time steps can cause rotation matrices to lose their 
orthonormal properties. Euler angles suffer from gimbal lock, where at certain orientations, two 
Euler angles describe rotation about the same axis, and a degree of freedom is lost. Oriented 
angles are inconvenient to accumulate (in particular one cannot simply add oriented angles to 
represent successive rotations) and require the evaluation of trigonometric functions. In this work, 
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we choose to use normalized quaternions, which require 4 floating point numbers to store, are easy 
to normalize, can be accumulated in a convenient manner, and avoid the need for (potentially 
expensive to evaluate) trigonometric functions. 

A normalized quaternion can be used to represent a finite rotation relative to a given initial 
reference frame, and is specified by 9 = {s,p} £ R , a combination of a scalar s and a vector 
p £ M 3 that satisfy the unit-norm constraint 

|| 0|| 2 = s 2 + p ■ p = 1 . ( 1 ) 


Quaternions can be combined via the operation of quaternion multiplication, whereby 9 3 = 6162 
is defined via 


S 3 


SlS 2 -P 1 P 2 

P3 


SlP 2 + S 2 P 1 + Pi x p 2 


( 2 ) 


with 9i = {sj, Pj}, i = 1,2,3. With this operation, normalized quaternions form a group with 
identity { 1 , 0 }; the inverse of a quaternion 9 = {s,p} is given by 9 1 = {s, — p}. 

In this work, we will use normalized quaternions to represent the orientation of a body in three 
space dimensions. Any finite rotation can be defined by its oriented angle, a vector 0, indicating 
a turn of 0 = || 0 || radians counterclockwise (i.e., using the right-hand convention) around an axis 
0 = 0/0. This rotation can be associated with the quaternion 


9$ = (cos( 0 / 2 ), sin( 0 / 2 ) 0 }, 


(3) 


i.e., p gives the axis of the rotation and the magnitude of p gives the angle of rotation; the inclusion 
of s and the normalization constraint is thus not strictly necessary [33] but is useful numerically. 
Note that 9 and —9 correspond to the same physical rotation/orientation [BSj. 

Performing a rotation on any three dimensional vector r in the reference frame gives a rotated 
vector r' = R(9)r, where the rotation matrix is 

pp T + sP + ^s 2 - 0 I . 

Here P is a cross-product 3x3 matrix such that Pr = p x r for any r, i.e., Pij = eikjPk , where 
e is the Levi-Civita symbol. Given two normalized quaternions 9\ and 9 2 , their rotation matrices 
satisfy the condition 


R(9) = 2 


i*(0i)i*(0 2 ) = R(0i9 2 ) 


(4) 
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that is, successive rotations can be accumulated by multiplying their associated quaternions. More 
precisely, if a rotation given by oriented angle <fi followed by a rotation -0 yields a total rotation £, 
then it holds that 6 ^ = 6^9^. 

Given an angular velocity u>, we can write the corresponding time derivative of orientation as 

6 = T'cj, (5) 


where th (6) is the 4x3 matrix 


^ = - 

2 


~P 

si -P 


( 6 ) 


The matrix SI/ has many properties that will be useful when we formulate equations of motion for 
bodies with orientation. First, it satisfies the property 


\ {~sp + sp) = 0, (7) 

which together with the relation 6 = indicates that the deterministic evolution (JH]) remains 
on the constraint 0- This property is used in Section [TT] to show that the Langevin equations 
presented in this work also preserve the constraint. Another useful relationship is the fact that 


de ■ tF 7 = 0 i.e. di (4! ik ) = 0, 


( 8 ) 


which is clear because the j-th row of \F has no entries that depend on the j-th component of 6 . 
Here and in the remainder of this paper we use Einstein’s repeated index summation convention, 
and denote dj = d/dOj. 

Describing the orientation of a body at several times t n requires choosing a single initial reference 
orientation associated with 9° = {1,0}, and recording the quaternion G n that describes the rotation 
from the reference orientation to the orientation at instant t n . Furthermore, if the body undergoes 
a rotation with constant angular velocity u from time t n to time t n+1 = t n + At, we have that 
9 n+1 = O^AtO 71 . This leads to a natural recipe for tracking orientation using the Rotate procedure 

m 

6 P+ 1 = Rotate(0 n , uAt) = 9 wAt G n . (9) 


In constructing numerical schemes in Section III it will be necessary to consider the second order 
expansion of this rotate procedure 

(a; ■ lj) At 2 


Rotate(0, uAt) =9 + H/cjAf — 


9 + 0 (At 3 ) , 


( 10 ) 


as shown in Appendix (A 1). 
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B. Rotational Brownian Motion 

For simplicity, we first consider a single rigid body that is free to rotate but with a reference 
point q, around which torques are measured, that is fixed in space. We let the orientation of this 
body (relative to some fixed reference frame) be denoted by the quaternion 6 (t) , and we suppose 
that the body is subjected to a torque r generated by a given conservative potential U{6). It can 
be shown (see Appendix |A 2[ ) the the torque generated by the potential is 

T = -V T dgU (11) 

In practice, it is not necessary to formulate U(G) and calculate — dU/dd to obtain the torque. 
Often is is much more convenient to calculate torque directly based on the geometries of the rigid 
bodies and the forces applied to them. We will see that © will be a convenient relation for 
formulating the constrained equations of motion. The schemes that we develop will be able to 
simulate the motion of rigid bodies without direct knowledge of U(0 ); they simply update the 
positions and orientations of the bodies based on the total forces and torques applied to each body. 


1. Overdamped Langevin Equation 


We introduce the 3x3 symmetric positive semidefinite (SPD) rotational mobility matrix 
which acts on torque to produce the resulting angular velocity, u> = M^ t t. Note 
that the mobility contains all the effects of hydrodynamics, including the shape of the body, the 
hydrodynamic interactions with other bodies or boundaries, far-field boundary conditions, etc. In 
this section we will assume this matrix is known, and discuss ways to obtain it explicitly in Section 
IV Using ([5]) and ©> we can write down a deterministic equation of motion for the rigid body, 


— =^M WT r = - d e U = - MdgU , 

where we have defined the 4x4 mobility matrix M = M ut-4 1:1 . 

It is now straight forward to formulate an Ito Langevin equation for the rotational diffusion of 
the body, 


UQ _ _ 1 _ 

— = - MdgU + y/2k B T M 2 W + {k B T)do-M, 


( 12 ) 


where W (t) is a collection of independent white noise processes. Here M 2 = with the 


i i 

“square root” of the mobility obeying the fluctuation-dissipation relation ( M 2 


1 

2 

LOT 
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M, 


UT , for example, it could be the Cholesky factor of Note that in (12) the covariance 


for the noise satisfies the fluctuation dissipation balance condition M 2 ^M 2 ^J = M. The i-th 

component of the stochastic drift term dg ■ M may be written in indicial notation as djMji(G). 


Using Ito’s formula, we can show that the overdamped dynamics (12) strictly preserves the 
constraint that 6 have unit norm, 


J ( e T e ) =0 Td ^ + ( k B T ) I :M = ( k B T) (q t (dg ■ Af) + I : Af) = (k B T) dg • (6 T M ) = 0, 


dG 


where we used ([7| and its consequence G 1 M = 0. Note that the stochastic drift term in (|12|) can 
be rewritten as (see Appendix [B|), 


dg ■ M = dg ■ = * ( d e M ur ) : - -Tr (M ur ) G, 


(13) 


where Tr denotes trace, and colon denotes double contraction; in index notation 
(VidgM^) : * T ). = 'S'ijdt (M UT )j k ^>ik and {Tr (M u}T )G} i = (M tJT ) jj di . We will see that 


this way of writing the drift is convenient when we consider numerical methods for integrating (12) 


in Section III Note that the stochastic drift term proportional to Tr (M^t) 0/4 can be seen in Eq. 
(36) in [26] to be related to enforcing the normalization constraint, and it will turn out we do not 
need to include it explicitly just as in [26] . 


In the special case of a free particle with unit mobility, M WT = I , (12) degenerates to the 


Stratonovich equation (see (26)) 


0 = (2k B T)i <S>o W. 


(14) 


Recall that the infinitesimal change in orientation is given by the infinitesimal rotation dcj) in the 
axes-angle representation, where the direction of the vector dtp is the axes around which the body 
is rotated by an angle dp. Also recall that the corresponding change in the quaternion is 

dG = V (G) dp, 

at least deterministically. Since the standard rules of calculus apply in the Stratonovich interpre¬ 


tation, (14) is equivalent to 


dip = (2 k B T) 2 dB 


(15) 


where B (t) is Brownian motion, formally W = dB/dt. This is a natural definition of isotropic 
rotational diffusion [53]. 
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We can verify that (12) has the correct noise covariance when M.^ T is not a multiple of the 


identity by considering the rotational mean square displacement at short times. Let us consider 
a set of orthonormal vectors u^(t) which are attached to the rigid body, and define a rotational 
displacement following Kraft et al. [1], 

(16) 


1 3 

A u (At) = - ^2 u i(°) x u i (At) ■ 


i=l 


A straightforward calculation relates this rotational displacement to the total angle of rotation 
relative to the the initial configuration, 


Au (At) = sin (<j>At) &At = 4>At + °( At2 ), 


(17) 


which shows that the magnitude of the rotational displacement is insensitive to the choice of the 


initial triad u,;(0). If the covariance of the noise in (12) is correct, it should hold that (c.f. Eqs. 
(1,2) in Ref. gj) 

\T\ iT 


1 ((Aii (At)) (Au (At)) J 

2k B T At —>0 


1 


At 


2k bT At^o 


Urn ( 0A ' 0At 1 =M„ 


At 


(18) 


which follows directly from (12). This shows that our equation has the same physical noise co¬ 


variance as the overdamped equation in Ref. g|, only written in a different representation. In our 
numerical tests, we will use ((Aii (t)) (A u(t)) t ) as a convenient definition of a rotational mean 
square displacement (RMSD) at time r; note that this RMSD is necessarily bounded and thus 
must reach a plateau at long times. 


2. Smoluchowski Equation 


A key property of the overdamped Langevin equation (12) is that it is time reversible with 
respect to the Gibbs-Boltzmann equilibrium distribution 


Peq ( 0 ) = ^ exp (-U (6) /k B T) 5 (e T e - 1 ) , 


(19) 


with Z a normalization constant. The overdamped equation (12) has the familiar structure of 
a generic Langevin equation (see Section I.A in Ref. my, however, a crucial difference is that 


(12) is an SDE on a manifold, namely, the unit 4-sphere, rather than an SDE in Eucledian space. 


A discussion of overdamped Langevin equations constrained on a manifold can be found in Ref. 
[33]. As explained there, for general curved manifolds one has to carefully construct the stochastic 
drift terms in order to ensure consistency with the desired equilibrium distribution. Note that the 
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original (true or physical) dynamics is unconstrained, and could, in principle, be described using 
a non-redundant parameterization of the rotation group such as Euler angles [lOj; the unit norm 


constraint implicit in (12) arises because it is mathematically simpler to embed the unit 4-sphere in 
a four dimensional Euclidean space than to parameterize it directly. The geometric matrix (6) 
plays the role of the projection operator in Ref. [34) . but unlike a projection operator, T' is not 
square and projects from the original (physical) three-dimensional space of angular velocity to the 
tangent space of the unit 4-sphere. 


To see that (19) is indeed the equilibrium distribution let us consider the case of a freely- 
rotating particle, U (0) = 0, which should correspond to uniform probability of all orientations. 
The uniform distribution over the space of orientations of a rigid body in three dimensions is the 
so-called Haar measure over the group SO( 3), and has been the subject of mathematical study 
[35) [36] . It is known that in the quaternion representation this Haar measure corresponds to a 
constant probability density over the surface of the unit 4-sphere, i.e., the Hausdorff measure on the 
unit 4-sphere [331E51I37]; generating random uniformly-distributed orientations amounts to simply 
generating a point uniformly sampled on the unit 4-sphere m • This uniform distribution over the 


unit quaternion sphere is captured in (19) by the term 6 (O' 0 — l), and the additional prefactor 
exp (—17 (6) /ksT) captures the standard Gibbs-Boltzmann weighting of the configurations based 
on their potential energy. 

Note that more generally, for a manifold E defined by the scalar constraint g(0) = 0, the 
Hausdorff measure da s on the the surface contains a metric factor relative to the Lebesque measure 
dO in unconstrained coordinates, as given by the co-area formula 


dax(0) = 6(g(0))\\Vg(0)\\ 2 d0. 

In our case, however, g(6) = 0 1 0 — 1 and ||V<?(0)|| 9 = ||0|| 2 = 1 is constant over the surface of 
the unit 4-sphere, and the metric factor can be absorbed into the normalization factor Z. The 
fact that no metric factors appear in the quaternion representation simplifies the equations; in 
other representations such as Euler angles or rotation angles the Gibbs-Boltzmann distribution is 
not uniform even in the absence of external potentials, and therefore “metric forces” need to be 
included in the stochastic drift terms to ensure the correct equilibrium distribution [MUH]. This 
subtle point has been missed in a number of prior works even though the concept of metric forces 
is well understood for rather general constrained Langevin equations [38] • 


To demonstrate that (19) is the equilibrium distribution (invariant measure) for (12), we examine 
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the Fokker-Planck equation (FPE) for the probability density P (6,t), 

d t p = di [m 13 [(d 3 u) p + ( k B r) djP ]}. 


( 20 ) 


When P is the Gibbs-Boltzmann distribution (19), we formally obtain 


[(djU) P e q + (k B T) d 3 P e q] ~ exp (—17 ( 0 ) / k B T) 5' ( 0 T Q - l) 9,. 

We can then use the fact that \I/ 7 6 = 0 to see that at thermodynamic equilibrium the thermo¬ 


dynamic driving force inside the square brackets in (20) vanishes, which implies that the Gibbs- 


Boltzmann distribution is an equilibrium distribution; using standard tools combined with reason¬ 


able assumptions on U(0), it can also be shown that (19) is the unique invariant measure |33] - Note 
that the calculation above is formal, but one can make a more precise argument by considering 
the backward Kolmogorov equation applied to E [/] for an arbitrary well behaved function / and 
expressing expectation values as integrals over the unit 4-sphere, similar to the approach taken in 
J. A similar calculation can be used to show that the generator of the Markov diffusion process 


(12) is self-adjoint with respect to a dot product weighted by the invariant measure (19), which 


proves that the overdamped dynamics is time reversible with respect to (19). 


We can compare the Eq. (20) with the FPE derived for rigid rods in Ref. [5]. A rigid rod 


can 


be parameterized with a unit 3-vector ip indicating the orientation of the rod. If we expand (20) 
and use the property ([8]), we can rewrite the FPE in the form 

d t p =di (* ifc (M WT ) fcI [(djU) P + ( k B T) djP}} 

=V ik di {(M WT ) fct (Vji 0 djU) P + ( k B T) y 3 id 3 p )}. 

This FPE has exactly the same form as the rotational part of Eq. (4.149) in [5], with the crucial 
difference that for rods th is the cross product matrix corresponding to the direction ip. We see 


that (20) is a natural generalization of the standard Smoluckowski equation for uniaxial bodies to 


biaxial bodies. 


C. Rotation-Translation Coupling 

In order to describe Brownian motion of freely suspended particles, it is necessary to also include 
translation in our model of rigid body motion. We first consider tracking both the location and 
orientation of a single rigid body. To do this, we introduce a variable q ( t ) for the Cartesian 
coordinates of a chosen tracking point fixed in the body frame. We assume that we are given 
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hydrodynamic information in the form of a known grand mobility matrix N ( q,6 ), which is the 
linear mapping from given force F and torque r (about q ) to the resulting velocity u = q and 
angular velocity cj, 


u 

= N 

F 

_ 

u 


T 



uF Mi UT 

MuF 


( 21 ) 


where M UT = M^ F is the translation-rotation coupling tensor, and M u p is the translation- 
translation mobility familiar from Brownian Dynamics of spherical particles. 

Let us suppose that the torque and force are generated from a conservative potential U(q,6). 
Then using the fact that q = u, along with ([5| and © we can write the overdamped Langevin 
equation including translation as the Ito SDE, 


v 


d'T '—/ --—-—-- 

= — = - Nd x U + ^2k B TN 2 W + (k B T) d x ■ N 
at 

= - (SATH t ) d x u + y/2k B T3N*W + ( k B T) d x • (E1VH t ) , 


( 22 ) 


where x = ( q,6) and v = yu,6J are composite vectors of the translational and rotational 
variables (and their velocities), and we have introduced the block matrix 


I 0 
0 T' 


(23) 


The “square root” of the mobility N 2 satisfies the fluctuation-dissipation relation N 2 ( N 2 ) = 


AT". A similar computation to that mentioned in Section IIB shows that (22) is time reversible 
with respect to the Gibbs-Boltzmann distribution [33], 

(24) 


Peq(q , 0) =^ _1 exp (-U ( q , 0) / k B T) 5 {0 T 6 - l) 


III. TEMPORAL INTEGRATORS 


In this section we introduce temporal integrators for the overdamped equations of motion of 
rigid bodies immersed in fluid, as formulated in Section[nJ We update the quaternion representation 
of orientation using the Rotate procedure ([9]) introduced in Section IIA, preserving the unit-norm 
constraint to numerical precision. The stochastic drift term in ( |22| ) is approximated in two ways, 
using a Fixman midpoint scheme and a Random Finite Difference (RFD) scheme, see Section I.B 
in Ref. [2D] for a comparison of the two approaches in the context of unconstrained overdamped 
Langevin equations. 
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A. Euler-Maruyama scheme 


For illustration purposes, we begin by considering a naive Euler-Maruyama (EM) scheme applied 


to an incorrect variant of (12), in which we do not carefully handle the stochastic drift term 


(fcgT) dg ■ M. In the EM scheme, we advance the configuration from time level n to time level 
n + 1 with the time step 


u n = - M"r n + ( 2 W n 


v At 

6 n+1 = Rotate (6 n , u: n At), 


(25) 


where a superscript denotes the point in time at which a particular quantity is evaluated, e.g. 
6 n « 0(nAt) and 1W” T = M^t (0 n ), and the Rotate procedure is defined by ([£]). Here W n is a 
collection of i.i.d. standard (i.e., mean zero and unit variance) Gaussian variates generated using a 
pseudo-random number generator. Here and henceforth, we have used to express the updates 


directly in terms of torque t(6). While the scheme (25) is not actually consistent with (12), it 


makes a natural starting point when discussing temporal integrators for (12). 


Note that because we use the Rotate procedure this update actually moves along the 
unit norm constraint of normalized quaternions, as opposed to stepping off of the constraint and 
then projecting back onto it |3lj. This is a natural way to update orientation multiplicatively 
while still being consistent with the additive Langevin equations formulated in Section [ITJ In the 
alternative approach followed in [26] one has to solve a quadratic equation (c.f. (15) in |26]) 

for a Lagrange multiplier to enforce the normalization constraint; while this avoids the use of 
trigonometric functions, it is difficult to make such methods second-order accurate. We can expand 


the Rotate procedure using the Taylor series (10) and truncate the result at first order in At, to 


obtain an expression for the leading order change in 6, 

e n+l -e n rn n (u n -u n 

n u n - At- 


- 0 ". 


At ~ ~ 8 

='F n (-M” T r n + (2 k B TMl T f* W 7 
(W n ) T M n W n 

- ( k B T) y J “ T e n + o (At). 

Note that the last term is equal in expectation to — ksT (Tr (1W” T ) /4) 6 n , which gives us the 


second term in the stochastic drift on the right hand side of (13). Therefore, when constructing 


temporal integrators that are actually consistent with the correct dynamics (12), we see that we 


only need to add terms in the orientational update 0 n+1 — 6 n that will generate the remaining 
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stochastic drift (doM WT ) : ty 1 Fortunately, adding this term to the orientation looks 

to first order like a Rotate procedure with angular velocity k B T : VF 7 )"■ With this in 

mind, we now construct first order weakly accurate temporal integrators for 


B. Midpoint Scheme 


The standard approach to handling the stochastic drift in overdamped Langevin equations is to 
use Fixman’s midpoint scheme [391140| . Henceforth we consider the full equations (22) including 


translation and rotational diffusion. To apply the Fixman method to (22) we rewrite (22) in a split 
Stratonovich-Ito form, 

dx dU 


f = - ^ + (2 k B T)t HAT o N~\ W, 


(26) 


where the terms after o are evaluated at the beginning of the time interval in the spirit of the 
Ito interpretation, while the terms before o are evaluated at the midpoint of the time interval, 
in the spirit of the Stratonovich interpretation. Here 1V~ 5 satisfies IV - 5 ^_/V~ 2 ^ = JV“ 1 ; the 

term N~ 2 W can be thought of as a “random force and torque” fH] and is equivalent in law to 

of the Appendix. 


N l N 2 W. We demonstrate that (|26l) is equivalent to (|12|) in section 


B 1 


Note that the Fixman scheme can be seen as a direct application of the Euler-Heun nn 
predictor-corrector method [32] to the split Ito-Stratonovich form (26). We also ensure that the 


scheme is weakly second-order accurate for the linearized Langevin equations (i.e., for additive 
noise) by following the predictor-corrector approach described in detail in Ref. [25], giving our 
midpoint predictor-corrector Fixman-like temporal integrator, 


v n = (u n ,u n y = N 


p,n+i _ n 


T 

At 

2 " 


+ 




VF”’ 1 


(27) 


6 F ’ n+ 2 =Rotate ( B n . — 


-u 


V 


P’ n + 2 = N 


2 

P,n+\ 


+ 


k B T 

At 


jV p,n+ 5 (VF”’ 1 + W n ’ 2 ) 


g n+l =q n + A tv P,n+\ 

6 n+1 = Rotate (6 n , Atu p ’ n+1 2 


We show that this scheme produces the correct stochastic drift in Appendix [Bj more precisely, the 


scheme (27) is a first-order weak integrator for (22). 
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The Fixman scheme requires the application of iV _ 5 ; or, equivalently, of A/" -1 ; this is compu¬ 
tationally expensive in cases when only N is easy to compute, and it is prohibitive in cases when 
only the application of N and N 2 to a vectors can be computed. In the examples we study here 
these matrices will be small and thus easy to compute using direct linear algebra, but this approach 
does not extend easily to suspensions of many rigid particles. In the next section, we show how to 
avoid using N~ 2 or iV” 1 by using a random finite difference (RFD) approach. 


It is important to observe that the Fixman scheme (27) is unaffected by the change of represen¬ 
tations of orientations. All that needs to be changed to use other representations of orientations is 
to simply change the Rotate procedure. This point has already been intuited in prior works, where 
the standard Fixman method has been used for non-spherical bodies, such as, for example, work on 
Brownian dynamics for rigid rods m- The analytical simplicity of the quaternion representation 
makes it straightforward for us to prove first order weak accuracy for the Fixman scheme in the 
general case (see Appendix [Bj , although the key idea is in fact to write the dynamics in the split 
Ito-Strato form 


C. Random Finite Difference scheme 


To avoid the computation of AT 2 or N 1 , we formulate a random finite difference (RFD) 
scheme H51 [29] by expanding the stochastic drift term into pieces (see Appendix [B]) , 

0 

{ds'&it) Ik su 

where i,n,m, and p represent any component of x, and s,t, and u represent indices that range 
over only the orientation components, i.e. components of 6. 


{d x -(SNS T )}.=E im ( d n N, 


mp ) i —'np 




(28) 


An Euler-Maruyama scheme such as (25) will, in expectation, produce the second term on the 

The remaining term E ( d x (N) : S 1 ) can be 


right-hand side of (28), as we saw in Section 


III A 


approximated in expectation using an RFD correction [15J [29] to the velocity as follows, 

. r p ' n 

Y =W (29) 


v = ( 

u n , u 

x = ( 

> «) 

n ) = 

1 N 


F 

T 


+ 


k B T 


2k ^ T - (N^Yw n +'^(N(x)-N n )W" 

At V / 0 


^n+1 


q- - =q n + Atu n 
0 n+1 = Rotate (6 n , Atu : n ), 
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where W is a collection of i.i.d. standard normal variates generated independently at each time 
step. Here 5 is a small parameter that should be chosen to minimize roundoff errors in the finite 
difference (29] . Observe that the RFD scheme only requires the application of TV n , Ala , and N(x), 
which can be a considerable advantage over the Fixman scheme in the case when TV is expensive 
to invert. In Appendix [B] we show that this RFD approach does in fact generate the correct drift 
terms, and the scheme (|29l) is a first-order weak integrator for (22). Observe that the RFD scheme 


is also invariant under changes of representation for the orientation of the body; all that is required 
is an appropriate Rotate procedure. 

D. Suspensions of Rigid Bodies 


The temporal integration schemes presented above straightforwardly generalize to suspensions 


of more than one rigid body. The overdamped Langevin equation (22) continues to hold, but now x 


collects the positions and orientations of all bodies, and v collects the linear and angular velocities 


of all bodies, and E is a block-diagonal matrix with one diagonal block (23) per body. We assume 


here that we can compute the grand mobility tensor TV for all of the bodies, which maps the forces 

and torques applied on the bodies to the resulting linear and angular velocities, and accounts for 

the hydrodynamic interactions among the bodies. 

The deterministic and stochastic terms are handled in a straightforward way; we accumulate 

deterministic velocities and angular velocities on each body using the grand mobility tensor, and the 

random velocities and angular velocities that the bodies experience are given by yj2ksT / Nt TV 2 W 
.INI’ 


where TV 2 ^_Z\7 2 J = TV. Note that the direct computation of TV 2 can be expensive in the multi- 
body setting; a generalization of the fluctuating immersed boundary method m or the fluctuating 
force coupling method [H] can, however, generate the stochastic forcing in essentially linear time 
by using a fluctuating hydrodynamic solver. 

We focus here on generalizing the Fixman and RFD approximations of the stochastic drift 
term. The grand mobility TV consists of blocks Nab which take forces and torques on body B and 
produce the resulting velocities and angular velocities on body A (which can be the same as body 
B). We consider now the stochastic drift for a given body A, denoting the position and orientation 
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of body A with x A = {q A i &a}, 
dxA 


drift = (k B T) ^ d Xg ■ (S A Nab^b) > 

B 

= ( k B T) [( d XB 3 A ) : (JVabS£) + S A d XB {N AB S t b )] 

B 

= ( k B T) (8 XA S A ) : (N AA 3 t a ) + (fc B T) £ H a {d Xg N AB ) : 3%. 


(30) 


B 


where the sums range over all bodies B, and we used the fact that d XB S A is nonzero only when 


A = B. The first term term on the right hand side of (30) is a local term that does not contain 


any many-body effects, and can therefore be approximated by using the Rotate procedure, as for a 


single body. The second term on the right hand side of (30) can be approximated using a random 


finite difference or Fixman approach in the same way as for a single body. This second term 
contains many-body interactions which are captured in the computation of IV - 3 in the Fixman 
approach, and in the RFD approach they are captured by randomly displacing all bodies together 
(rather than one by one). 


IV. DIFFUSION ALONG A NO-SLIP BOUNDARY 

In this section we apply the Fixman and RFD temporal integrators to several model examples 
of a single rigid body immersed in a viscous fluid. Since we want to focus on examples with 
configuration-dependent mobilities, we examine rigid bodies confined to be in the vicinity of a 
no-slip boundary. Specifically, we simulate the diffusive motion of a tetramer of colloidal spheres 
(Section IVC), an asymmetric sphere (Section |IVD|, and a colloidal boomerang (Section |IVE[ ), in 
the presence of gravity and a no-slip wall located at the plane z = 0. For the tetramer we validate 
our new methods by comparing to the FIB method m with stiff springs used to keep the tetramer 
nearly rigid. The Python codes used to produce the results reported here are available as open 
source at https://github.com/stochasticHydroTools/RotationalDiffusion, 

There are two main types of quantities that we examine in these simulations, the first static and 
the second dynamic. The first type of quantities are various moments of the equilibrium distribution 
for the position and orientation of the rigid bodies, which we compare to moments of the expected 
Gibbs-Boltzmann distribution. The second type of quantities we study are components of the mean 
square displacement of the rigid bodies, as we now explain in more detail. 
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A. Mean Square Displacement 


D(t;x) =(A X(t-x) (AX(r; x)) T ) = 


(t;x), 


(31) 


We define the total mean square displacement (MSD) at time r as the outer product 

D t D 
Do D 

where A X(t-,x) = ( Aq(r-,x ), A u(t-,x)), with position increment Aq(r) = g(r) — g(0) and orien¬ 
tation increment Att(r) as defined in 


The average in (31) is taken over trajectories started at 


x = (q (0), 9 (0)). The short-time diffusion coefficient is given by the mobility in agreement with 
the Stokes-Einstein (SE) relation 

1 D(t;x ) 


lim 

2k bT t —»o r 


= N(x), 


(32) 


where the grand mobility tensor N is the block matrix (21). Our overdamped Langevin equation 


(22) is consistent with the SE relation (32), specifically, using © for the rotational component 
D r , we obtain 

(d T d^) (d r 0r) 

(<M^) (trtf) 

where d T = A q(r; x) is the translational displacement and </> r is the angle of rotation over the 


D(t-,x) = 


+ 0 ( 7 - 5 ), 


(33) 


short time interval r. The SE formula (32) follows directly from (33) the noise term in (22) and 


(18). 


We can further define the equilibrium MSD via the ergodic average 

D{t) ={D (t; x)), 


(34) 


where the average is taken over x distributed according to the Gibbs-Boltzmann distribution (24). 


In practice, we estimate D(t) from our simulations by taking a time average over one long trajectory 


(using the ergodic property) with the initial condition distributed according to (24), as generated 


using an accept/reject Monte Carlo method. We estimate error bars by running an ensemble of 
statistically independent trajectories. 


The Stokes-Einstein relation (32) gives the short-time mean square displacement, which can be 


used to define a short-time translational diffusion tensor 

Xst = l lim n ( ' Dt< ' T,X ^ = k B T ( M uF (*)). 

2 r—> 0 r 

In general, it is much harder to characterize the long-time diffusion coefficient 

x„ = i l.m 

2 t —>00 r 


(35) 


(36) 
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in the presence of confinement, even for a single body. The only simple case is when the MSD is 
strictly linear with time so that the long and short time diffusion coefficients are equal and one can 


just average the mobility over the Gibbs-Boltzmann distribution in order to obtain Xit using (35). 

Observe that the long-time diffusion coefficient is independent of the choice of point to track on 
the body, i.e., the choice of the point around which torques are expressed, however, the short-time 
one does depend on the choice. Our goal will therefore be to identify a tracking point that makes 
the MSD as close to linear as possible, so that the short-time diffusion coefficient provides a good 
estimate of the long-time one. If this can be accomplished, then the long-time diffusion coefficient 


can be estimated from the purely equilibrium calculation (35), without requiring us to simulate 


long trajectories and use (36). 


1. Choice of tracking point 

The choice of the origin around which torques are expressed, which is the point on the body 
whose position we track when computing the translational MSD, strongly affects the short-time 
MSD. Given two fixed tracking points q x and q 2 , it is straightforward to derive the following 
relationship between the blocks of IV 1 calculated using origin q ll and N 2 calculated with origin 

q 2 EEL 

Mi T =M\ yr (37) 

M'IjF =M) jjF + X T"12 

M 2 uF =M 1 uF - n 2 x (Ml T x r 12 ) + (M^ f ) T x r l2 - r l2 x M^ F , 

where r± 2 = q 2 ~ Qi- Cross-products between vectors and tensors are defined in Eqs. (4,5) in (43j, 
with A x b corresponding to taking cross products between rows of A and 6, in index notation, 

(A x b) {j = (Ai ■ x b) j = e jk iA ik bi , (38) 

where e is the Levi-Civita tensor, and similarly, 

(b x = (b x A :> j) i = e ik ib k Aij. (39) 

In general, the cross-coupling (translation-rotation) mobility tensors M^f = are not 

symmetric. However, it can be shown that for any body shape, there exists a unique point in the 
body called the center of diffusion or center of mobility (CoM), such that, when that point is taken 
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as the origin, the coupling tensors are symmetric, M^ F = M^f = M UT = MJ it . The location 
of the CoM can be found by solving for i j the linear system 


e ikl 


e jkl {Mu> T )ik 


CoM 

r l 




(40) 


where the mobilities are evaluated at an arbitrary origin, and r CoM goes from the origin to the 
CoM. It is very important to note that in the presence of confinement the location of the CoM is 
not fixed relative to the body but changes with the position and orientation of the body relative to 
the boundaries. Therefore, one should consider the CoM computed in the absence of confinement 
only as an approximate CoM. For some bodies of sufficient symmetry, there exists a point called the 
center of hydrodynamic stress (CoH) [44] . where the cross-coupling vanishes, M^f = M UT = 0; 
note that if a CoH exists it is also the CoM. A CoH always exists in two dimensions, however, for 
general skew bodies in three dimensions a CoH does not exist |43|. In two dimensions the CoH 
is the origin for which a torque out of the plane does not induce any translational motion in the 
plane, and its position can be found from m 


r CoH = (-MupJMur, MvpjMur). 


(41) 


Experimental investigations have lead to the suggestion that for planar particles confined to perform 


essentially quasi two-dimensional diffusion, the point (41) should be tracked mm- 


2. Free isotropic diffusion 


For a freely-diffusing rigid body in an unbounded fluid in the absence of any external forces and 
torques, all orientations are equally likely. It is well-known that in the oriented angle representation 
the Haar measure over the rotation group corresponds to 0 uniformly distributed over the unit 
3-sphere, and a probability density P{4>) = (2/7r)sin 2 (0/2) for the angle of rotation (see (14) in 
[36] ). Combined with © this shows that for free isotropic rotational diffusion the asymptotic 
long-time value of the rotational MSD is finite, 


lim D t (t) 

T—>• OO 


= lim ((Am(t)) (A u(t)) T ) = 

T—¥ OO \ / 



-I 

6 


0.167 1, 


2 
7r 



sin 2 (0/2) sin 2 (0) d<f> 


(42) 


independent of the shape of the body. 

When the CoM is used as the tracking point, the translational MSD for free isotropic diffusion 
is strictly linear in time, 


D? oM (t) = ( 2 X t) • /, 


( 43 ) 
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where the average (short- and long-time) diffusion coefficient is \ = ( k B T) Tr(M u p) /d, with d 
being the dimensionality and the orientation of the body used to evaluate M u p being arbitrary. 
See discussion around Eq. (46) in [20] for more details, and note that the mathematical reason 
behind this result is the fact that for a symmetric cross-coupling mobility M UT , the translational 
component of the thermal drift term ( k B T ) dg • (A4" llT \I' i ) vanishes identically for a freely-diffusing 
body, see (31,32) in [213] . 


3. Confined Diffusion 

Our primary focus in this work is diffusion of rigid bodies in the presence of confinement and 
gravity, specifically, we consider diffusion in the vicinity of a no-slip boundary [35]. In a number 
of experiments, the Brownian particles being tracked are substantially denser than the solvent 
and thus sediment close to the bottom microscope slide due to gravity mm, or, the particles 
are confined in a narrow slit channel mm- In both cases the boundaries strongly modify the 
hydrodynamic response. Notably, the CoM will depend on the position of the body relative to the 
boundary, and for non-skew particles, there may not be a CoH in the presence of a boundary even 
if there is one in an unbounded domain. Note that in the presence of gravity the typical height of 
a rigid body above a plane wall is on the order of the gravitational height h g k B T/ ( m e g ), where 
m e is the excess mass of the particle relative to the solvent, and g is the gravitational acceleration. 
The value of the gravitational height varies widely in experiments from tens of nanometers to tens 
of micrometers, depending on the size and density of the colloidal particles. 

In the numerical studies that follow we examine the MSD of isolated rigid particles sedimented 
near a wall in the presence of gravity. We orient our coordinate system so that the x and y axes 
are parallel to the wall and the 2 axis is perpendicular to the wall located at z = 0. In experiments 
based on confocal or optical microscopy, only the motion of the particle parallel to the wall can be 
observed and measured, in particular, what is measured is the parallel mean square displacement 

D\\ (t) =D xx (t) + D yy (r). 

In our simulations, we apply no forces in the x and y directions, and at large times we expect that 
-D|| (r) will grow linearly with slope proportional to the long-time quasi-two-dimensional diffusion 
coefficient X 2 D which can be measured from simulations or experiments, 


-Dm (t) ~ at long times. 
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In general, we do not expect .Dm( r) to be strictly linear in time. However, if it is, then the long¬ 
time diffusion coefficient is the same as the short-time diffusion coefficient, and can be obtained 
by averaging the parallel translational mobility over the equilibrium Gibbs-Boltzmann distribution 

©, 

X2D = k B T <M,|) = k B T { M Bx f x ) = k B T ( Mp y p y ) . (44) 


If a CoH exists, and is independent of the configuration, then translational and rotational motion 


will decouple and (44) will be exact. In general, however, a CoH does not exist in three dimen¬ 
sions even in the absence of confinement. However, as argued in Refs. mi, if the diffusion is 
strongly confined to be effectively two-dimensional, either because of strong gravity or because of 
the presence of two tightly-spaced confining walls, an approximate CoH should exist and therefore 
the parallel MSD will be approximately linear in time. We will examine these claims numerically 
in Section IIV El 

We also investigate the perpendicular mean square displacement for the height above the plane 
wall, 


D_l(t) =D zz (t ), 


which we expect to reach a finite asymptotic value at large times due to the presence of gravity, 

lim D ± (t) ={{zi - z 2 ) 2 ) (z Z2 y (45) 

where z\ and z 2 are the heights of the tracking point for two configurations sampled uniformly and 


randomly from the Gibbs-Boltzmann distribution (24). A good generalization of the concept of a 
gravitational height for nonspherical particles is 


h g = [ \ lim D_ l (r) 


h 


f k B T 

\rn e g 


(46) 


where f„ is a geometric factor that is hard to compute analytically for a general body but can 


be computed using a Monte Carlo algorithm from (45); the factor of 1/2 is chosen so that h g = 
k B T / ( m e g ) for a point particle. Note, however, that f g depends on the choice of the tracking 
point, and should therefore be associated with a particular fixed point on the body. Also note 
that h g measures the relative displacement of the particle in the vertical direction rather than the 
distance to the plane; it may therefore be more appropriate to think of it as gravitational thickness 
rather than height. 
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For rotation, we examine the diagonal components of the R.MSD D r (r ), which must reach a 
finite asymptotic limit at large times since the rotational displacements are bounded, 


lim D v (r) = ((A u (0i, 0 2 )) (A« (01, 02)) T >(0 1 ,e 2 ), 


(47) 


where Au ( 61 , 62 ) is the rotational displacement (16) between two random orientations 0i and 0 2 


sampled uniformly from the Gibbs-Boltzmann distribution (24). 


In order to compute averages over the Gibbs-Boltzmann distribution, which includes the effects 
of gravity and steric repulsion from the wall, we use a Monte Carlo method to generate random 


samples distributed according to (24). The simplest way to do this is an accept-reject method in 


which we first generate a random position q with uniform height and a uniform random orientation 
0 of the body, and accept the random configuration x with probability exp (—U(x)/kBT). Note 
that in principle the Gibbs-Boltzmann distribution is unbounded in the z > 0 direction, and cannot 
be captured exactly by such an accept-reject method with height distributed in a finite interval. 
However, since the probability decays monotonically in the tail as ~ exp (—m^gh/ksT) where h 
is the height of the tracking point of the body, we can adjust the upper bound of the uniformly 
distributed height empirically to ensure that we are only neglecting an insignificant portion of the 
distribution. One can avoid this bias by using an exponential distribution as a proposal density 
in the accept/reject process, or by using Markov Chain Monte Carlo to generate samples from 
the Gibbs-Boltzmann distribution; we have found this to produce indistinguishable results for our 


purposes, while being significantly slower. We estimate asymptotic values of the MSD from (45) 


and (45) by using the Monte Carlo method to generate a large number of pairs of samples from 


the equilibrium distribution, calculating the mean square displacement between each pair, and 
averaging this value over all of the pairs. 

When using discrete time steps with stochastic forcing, it is possible for unphysical configu¬ 
rations to occur; this leads to a finite-time breakdown of explicit integrators mm such as our 
Fixrnan and RFD schemes. Specifically, in our numerical tests, it is possible for the stochastic 
terms to “kick” some part of the body through the wall; this invalidates the hydrodynamic calcula¬ 
tions used to compute the mobility, or makes the mobility not positive-semidefinite. To handle this 
possibility in our simulations, after each configurational update (including the predictor step to the 
midpoint in the Fixrnan scheme), we check whether any part of the body overlaps the wall, and if 
the new configuration is not valid, we start again at state x n and repeat the time step, drawing 
new random numbers. This procedure is repeated until a valid new state is found (note that it 
is possible for multiple rejections to occur in one time interval). Because this rejection of invalid 
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Figure 1: Rigid multiblob models of three types of particles studied in this work. The blobs are shown as red 
spheres of radius equal to the blob hydrodynamic radius, and the no-slip bottom wall is shown as a green 
plane. {Left) A spherical colloidal “surfer” that has a much denser metallic cube of hematite embedded in 
it, taken from the work of Palacci et al. (8]. In our computer simulations we model this as an icosahedron 
of rigidly-connected blobs, one of which (indicated by a blue sphere) models the dense hematite and holds 
all of the mass of the particle. {Middle) A tetramer formed by connecting four colloidal particles using 
DNA bonding into a tetrahedron, as in the work of Kraft et al. T]. The multiblob model has four blobs 
rigidly placed at the vertices of a tetrahedron. {Right) A right-angle boomerang colloid manufactured using 
lithography and studied in a slit channel formed by two microscope slides by Chakrabarty et al. P] , modeled 
here using a 15-blob approximation. 

states changes the dynamics, and therefore the statistics of the system, we ensure that the number 
of rejections is very low compared to the total number of steps taken. In the results presented 
in this section, the rejection rate (number of rejections divided by number of attempted steps) is 
never greater than 5 x 10~ 5 ; in most cases it is zero. 

B. Rigid multiblob models for hydrodynamic calculations 

For the purposes of hydrodynamic calculations, we discretize rigid bodies by constructing them 
out of multiple rigidly-connected spherical “blobs” of hydrodynamic radius a. These blobs can 
be thought of as hydrodynamically minimally-resolved spheres forming a rigid conglomerate that 
approximates the hydrodynamics of the actual rigid object being studied. Examples of such “multi¬ 
blob” |38j models of several types of rigid bodies studied in recent experiments are given in Fig. [lj 
As Fig. [2] illustrates for a rigid sphere, the hydrodynamic fidelity of rigid multiblob [45] models can 
be refined by increasing the number of blobs (and decreasing their hydrodynamic radius a accord¬ 
ingly); of course, increasing the resolution comes at a significant increase in the computational cost 
of the method. Similar “bead” or “raspberry” models appear in a number of studies of hydrodynam¬ 
ics of particle suspensions mma HEi-isi], with the blobs or beads being either connected rigidly as 
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Figure 2: Rigid multiblob or “raspberry” models of a rigid sphere, containing 12 blobs (left), 42 blobs (center), 
or 162 blobs (right) placed on the surface of a sphere using a recursive triangulation procedure starting from 
an icosahedron (left-most panel). The radius of the red spheres is equal to the effective hydrodynamic radius 
of a blob; the hydrodynamic radius of the resulting rigid multiblob sphere is computed numerically |48L I57| 
and is generally larger than the geometric radius of the sphere used for the recursive triangulation. 

we do here, or connected via stiff springs; in some models the fluid or particle inertia is included 
also. Since in this work we focus on the long-time diffusive dynamics it is crucial to use rigid 
rather than stiff springs, and to eliminate inertia in the spirit of the overdamped approximation, 
in order to allow for a sufficiently large time step to reach physical time scales of interest (seconds 
to minutes in actual experiments). 

After discretizing a rigid body using n blobs, we write down a system of equations that constrain 
the blobs to move rigidly. These intuitive equations are written in a large number of prior works 
UMaiMM EBl but we refer to j52j for a clear yet detailed exposition; the authors also provide 
associated computer codes (not used in this work) in the supplementary material. Letting A = 
{Ai,...,A n } be a vector of forces (Lagrange multipliers) that act on each blob to enforce the 
rigidity of the body, we have the linear system for A, u and u, 

^2(M B ) ij X j =u + uj x (n- q), Vi (48) 

3 

I> =*• 

i 

^(r*j - q) X A i =r, 
i 

where u is the velocity of the tracking point q , u is the angular velocity of the body around q, F 
is the total force applied on the body, r is the total torque applied to the body about point q, and 
r ? ; is the position of blob i. 
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Here the blob-blob translational mobility Mb describes the hydrodynamic relations between 
the blobs, accounting for the influence of the boundaries. The d x d block of the translational 
mobility ( Mb ) ij - computes the velocity of blob i given forces on blob j, neglecting the presence 
of the other blobs in a pairwise approximation. In the presence of a single wall, an analytic 
approximation to ( Mb)^ is given by Swan and Brady [35j|, as a generalization of the Rotne-Prager 
(RP) tensor [53] to account for the no-slip boundary using Blake’s image construction. In this 
work we utilize the translation-translation part of the Rotne-Prager-Blake mobility given by Eqs. 
(Bl) and (C2) in [55] to compute Mb , ignoring the higher order torque and stresslet terms in the 
spirit of the minimally-resolved blob model m- The self-mobilities for a single blob are given in 
(Dl) in Appendix [Dj Note that in a suitable limit of infinitely many blobs of appropriate radius, 
solving ( |48[ ) computes the exact grand mobility for the rigid body (or a collection of bodies), even 
though only a low-order RP-like approximation is used for Mb [571158]. To see this, note that blob 
methods can be considered as a discretization of a regularized first-kind integral equation [58] for 
the Stokes flow around the rigid bodies. In a recent experimental and computational study [60], 
the mobility of a rigid rod tethered to a hard wall was computed using the multiblob approach as 
we do here, however, in that work the wall was created from (many) blobs instead of using the 
known generalization of the RP tensor to a single-wall geometry [35] as we do here. 


The solution of the linear system (48) defines a linear mapping from applied force and torque 
to body velocity and angular velocity, and thus gives the grand mobility N (for explicit formulas, 
see [52]). Observe that generalizing the system (48) to a collection of rigid bodies is trivial. In the 


examples considered in the work, the number of blobs is small and the system (48) can easily be 
solved by computing the Schur complement [52] and inverting it directly with dense linear algebra. 


The use of dense linear algebra allows us to focus our attention on the temporal integrators for the 
overdamped dynamics and not on linear algebra or hydrodynamics issues. In principle, our tem¬ 
poral integrators can be used with a variety of methods for computing the hydrodynamic mobility 
of suspensions of rigid bodies, for example, boundary-integral or boundary-element methods can 
be used to compute the (action of the) grand mobility with higher accuracy. 

We compute the square root AT 5 by performing a dense Cholesky factorization on TV. It is 


important to note that if A 1b is SPD, the grand mobility N computed by solving (48) is also 
SPD. Note that the Swan-Brady approximation to Mb 


used here is based on the Rotne- 
Prager tensor and is thus only guaranteed to be positive definite when the blobs do not overlap 
each other or the wall, i.e., when no two blobs are closer than a distance 2a and the distance 
of all blob centroids to the wall is greater than a. It is possible to generalize the Rotne-Prager- 
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Yamakawa tensor to confined systems [61], thus guaranteeing an SPD Mb even when blobs overlap 
each other or the wall (but of course their centroids must remain above the wall), but we know of 
no published explicit formula that accomplishes this even for the case of a single wall. Fortunately, 
for our model of boomerang-shaped particles we numerically observe an SPD mobility when the 
blobs do not overlap the wall even though blobs overlap each other. 


C. Colloidal tetramer: Tetrahedron 


In this section we study a tetramer formed by rigidly connecting four colloidal spheres at the 
vertices of a tetrahedron EH], diffusing near a single no-slip boundary. The tetrahedron is dis¬ 
cretized in a minimally-resolved way using 4 blobs, one at each vertex of a regular tetrahedron, 
as illustrated in the pictured in the left panel of Fig. [l| In some arbitrary units, each blob is a 
distance d = 2 away from all of the others and has a hydrodynamic radius of a = 0.5; this some¬ 
what arbitrary choice makes the tetrahedron hydrodynamically sufficiently different from a sphere 
to require resolving the orientation of the body as well. 

To avoid symmetries and make the test more general, we assume each of the four spheres to 
have a different density; the gravitational forces on the vertices in the negative 2 direction are set 
to F\ = 0.15, F '2 = 0.1, F$ = 0.3 and F 4 = 0.05 in units of ksT/a. To prevent the tetrahedron 


from passing through the wall, we include a repulsive potential (49) between each of the blobs and 
the wall, based on an ad-hoc combination of a Yukawa and a hard-sphere-like divergent potential, 

h — a' 


^wall(^’ a ) 


eo 

h — a 


exp 


(49) 


where h is the height of the center of the blob above the wall, e = 20 ksT is the repulsion strength, 
and b = 0.5a to be the Debye length (these values are selected somewhat arbitrarily). The total 
force and torque on the rigid tetramer is the sum of the forces and torques on the individual blobs. 


The above choice of parameters gives the center of the tetrahedron a gravitational height (46) of 
h 0 ~ 1.75a. 


For comparison and validation, we also construct an approximation to the freely-moving rigid 
tetrahedron using four blobs connected by stiff springs, and then employ the FIB method [15] to 
simulate the diffusive motion of the almost-rigid tetramer; the same parameters are used in both 
simulations. The FIB simulation was performed in a domain of 64 x 64 x 64 grid cells of width 
Ax = 0.796a using the 4-point Peskin kernel, which ensures that the effective hydrodynamic radius 
of the blobs is a DS!- The boundaries on the top and bottom of the domain are both no-slip 
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walls, which differs from the domain for the rigid body simulations, but since the center of the 
tetrahedron almost never goes past 20% of the channel width (see figure [3]) , the effect of the top 
wall is relatively minor. The spring stiffness was set to k = 200 ksT/a to keep the deformations of 
the tetrahedron small; this imposes a stringent limit on the time step size At. 


1. Equilibrium Distribution 


In this section we examine the equilibrium distribution for the colloidal tetramer. We use a 
Monte Carlo method to generate the marginal Gibbs-Boltzmann distribution for the height of the 
geometric center of the tetrahedron, and compare to our numerical results. We see in Fig. [3] 


that the Fixman (27) and RFD schemes (29) are in good agreement with the Gibbs-Boltzmann 


distribution. The Euler-Maruyama scheme (25), with the obvious additions to include translation, 


however, neglects parts of the stochastic drift and generates an equilibrium distribution which has 
clear errors that do not vanish as the time step size is refined (not shown). 


2. Mean Square Displacement 


In this section, we examine the translational mean square displacement of the tetrahedron. In 
the left panel of Fig. [4] we examine the effect of the choice of tracking point on the parallel mean 
square displacement by comparing Du (r) when tracking the geometric center of the tetrahedron, 


versus tracking one of the four vertices. In both cases (32) gives the initial slope of the MSD as it 
must, and these slopes are clearly different. Since at long times the slopes of the parallel MSD is 
independent of the choice of tracking point, the MSD cannot be linear at all times for both choices 
of tracking point. Indeed, the results in Fig. [4] show that the parallel MSD is linear to within 
statistical and numerical truncation errors only when the geometric center is tracked. By contrast, 
the rotational MSD is insensitive to the choice of tracking point, as seen in the right panel of Fig. 

II 


We note that far from the wall, torques applied about the center of the tetrahedron generate 
no translation, indicating that in the absence of confinement the geometric center is both the 
CoM and the exact three dimensional CoH (which does not exist for general rigid bodies). In the 
presence of the boundary, this is not strictly the case, but we nonetheless observe in Fig. [4] that the 
average parallel mobility evaluated using the center of the tetrahedron as an origin gives a good 
approximation to the long time quasi two-dimensional diffusion coefficient \ 2 D • This is perhaps 
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Figure 3: Equilibrium distribution for the height of the geometric center of the tetramer colloid pictured in 
the left panel of Fig. [l] The results obtained by using the FIB method and the Fixman (Section III B) and 
RFD (Section III C|) integrators agree with the Gibbs-Boltzmann distribution generated using Monte Carlo 


sampling. The results obtained by using the inconsistent EM scheme (25) demonstrate that neglecting the 
stochastic drift term yields an incorrect equilibrium distribution. This plot is based on 16 runs of 3 • 10 5 
time steps with a small time step size of At w 0.0653 (6nr]a 3 /ksT ); no rejections were needed for this small 
time step in any of the runs. 


not surprising due to the high symmetry of a tetrahedron, as the geometric center is the “obvious” 
point to track. 

In Figure [5] we compare results for the MSD obtained using the overdamped rigid-body inte¬ 


grators from Section III to results obtained using the FIB method and stiff springs. We examine 
the mean square translational displacement parallel and perpendicular to the wall, as well as the 
rotational MSD, and find that the behavior of the tetrahedron is the same for both the stiff and 
rigid simulations; this provides a validation of our rigid-body methods and our codes. However, 
due to the presence of the stiff springs, using the FIB method to simulate a rigid body requires a 
time step size that is 32 times shorter. Due to the small time step size required for the tetrahedron 
constructed using rigid springs, and the high cost of numerically solving a Stokes problem each 
time step, it is computationally impractical to study the long time diffusion coefficient using the 
FIB method. The time step size for the rigid-body method could in principle be even larger and 















31 



0.25 


0.20 


: 0.15 


i 0.10 


0.05 


0 . 00 , 




"‘in, 


'***, 






RFD rot MSD vertex 
RFD rot MSD center 
asymp rot MSD 


100 


200 300 

Time 


400 


500 


Figure 4: Mean square displacement for a colloidal tetramer sedimented near a bottom wall. The data 
for both figures is generated from 4 independent runs of 3 • 10 5 time steps with a time step size of At w 
0.136 (6Trr]a 3 /ksT ); the highest rejection rate was 2.33 x 10~ 5 (a total of 7 rejections). The MSDs for each 
tracking point are calculated from the same trajectories. (Left) Comparison of parallel translational MSD 
D || (r) when tracking the geometric center (green squares), versus when tracking the fourth vertex (blue 


squares). We see that at short times the numerical slope agrees with (32), shown as a red dashed line of slope 
l.M(kBT/6nr)a) ~ 2.8 • 10~ 2 for the geometric center, and as a red solid line of slope 2 .21(kBT/6irria) « 
4.7 • 10 -2 for the vertex. (Right) Comparison of the parallel ((x — x) or (y — y)) component of the rotational 
mean square displacement using the two choices of tracking point. The asymptotic rotational MSD predicted 


by (47) is shown as a dashed line. 


still resolve the dynamics of the body, but it is limited by the stiff potential used to repel the 
particle from the wall; we keep At sufficiently small to strictly control the number of rejections 
of unphysical states where a blob gets too close to or passes through the wall. In Section [V] we 
discuss some ideas that may allow for the use of larger time step sizes even in the presence of steep 
repulsive forces. 


D. Asymmetric sphere: Icosahedron 

In this section we examine the diffusive motion of a rigid sphere whose center of mass is displaced 
away from the geometric center, in the presence of gravity and a bottom wall (no-slip boundary). 
This models recently manufactured colloidal “surfers” that become active when the particles sedi¬ 
ment to a microscope slide [8]; here we consider a passive particle in the absence of chemical driving 
forces. Diffusive and rotational dynamics of a symmetric patterned (Janus) sphere near a boundary 
has been studied experimentally by Anthony et al. [2j, and can be described well by theoretical 
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Figure 5: Comparison of the mean square displacement for a colloidal tetramer sedimented near a bottom 
wall, obtained by treating the body as rigid using the RFD and Fixman methods developed here (see caption 
of Fig. [4|, versus keeping it nearly rigid with stiff springs and using the FIB method m- For the FIB runs 
we used 32 simulations of 10 5 time steps each, with a time step size 32 times smaller than in the rigid-body 
simulations. (Left) Parallel translational MSD when tracking the geometric center of the tetrahedron. The 
inset focuses on the short time diffusion, and shows a slight hydrodynamic difference between the rigid and 
semi-rigid models that is due to the different methods used to handle the hydrodynamics, as well as the slight 
flexibility of the tetrahedron in the FIB simulations. (Right) The parallel ((x — x) or (y — y)) component of 


the rotational MSD (18). 


approximations for the mobility of a rigid sphere near a planar wall Tij . 

We construct a hydrodynamic model of an asymmetric rigid sphere of radius aj by rigidly 
constraining 12 blobs at the vertices of an icosahedron; a similar blob model of a sphere was used 
in Ref. [ 38] but was based on (stiff) penalty springs rather than rigid-body constraints. Note 
that more accurate results can be obtained by using more blobs to construct the spherical shell 
m ■ Each blob has a hydrodynamics radius of a = 0.175 and is located a distance 2.5a from the 
center of the icosahedron, so that the minimal distance between two blobs is about 2.63a. These 
parameters are chosen so that the icosahedron is hydrodynamically nearly rotationally invariant, 
and has an effective translational hydrodynamic radius (computed numerically) in bulk (i.e., far 
from the wall) of 


1 

67r?7 (Tr ( M uF ) /3) 


2.86a = 0.5, 


in some arbitrary units. A gravitational force of F = 0.5 = 1.25 ksT/ai is applied to one of the 12 
blobs, which represents the dense hematite cube embedded in the nearly spherical colloidal surfers 
of Palacci et al. [8]. Gravity therefore generates a torque around the center of the sphere and 
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causes the icosahedron to prefer orientations where the heavy blob is facing down. A short-ranged 


repulsive force U (h\ aj) given by (49) is added to keep the icosahedron from overlapping the wall, 
where now h is the distance from the center of the icosahedron to the wall, the repulsion strength is 
e = 20 ksT, and the Debye length is (arbitrarily) set to b = aj. This choice of parameters gives the 
center of the icosahedron a gravitational height ( |46[ ) of h g ~ 0.96a/. Note that in this example the 
icosahedron is considered to be a hydrodynamic approximation of a physical sphere and therefore 
the repulsive force acts on the center of the sphere (thus not generating any torque), rather than 
acting on each of the 12 blobs individually (which would generate some small spurious torque). 


1. Equilibrium Distribution 

We first investigate the equilbrium distribution P eq (q , 0 ), examining the marginal distributions 
of height h and orientation angle 6 , which is the angle between the 2 axis and the vector connecting 
the center of the icosahedron to the blob to which we apply the gravitational force. In this simple 
example, we can compute the marginals of the equilibrium Gibbs-Boltzmann distribution analyti¬ 
cally for both h and 6, and they are compared to numerical results in Fig. [6j We see that the RFD 
and Fixrnan schemes agree with each other and with theory. Due to the nonuniform gravitational 
forcing on the icosahedron, it prefers orientations with 6 closer to zero, but the thermal fluctuations 
causes it to explore all orientations. 

2. Mean Square Displacement 


To validate how well our scheme captures the dynamics of the system, we examine the mean 
square displacement of the geometric center of the icosahedron. In Fig. [7] we compare our results 
to the mean square displacement of an actual hard sphere with hydrodynamic radius a = 0.5. We 
apply torques and forces to the sphere that are identical to those applied to the icosahedron, but for 
the hydrodynamic mobility of the sphere we use the most accurate theoretical expressions available 
in the literature, see QD2|D3|D4D in Appendix [Dj instead of relying on the blob approximation to a 
sphere (Dl), even though in this specific case is sufficiently accurate. This tests allows us to 
both evaluate our temporal integration method, as well as to examine how well the 12-bead model 
approximates a single spherical particle. 

The results shown in Fig. [7] demonstrate that the dynamics of the icosahedral rigid multiblob is 
essentially identical to that of an actual sphere. Note that for a sphere the mobility does not depend 
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Figure 6: Equilibrium distribution for a rigid icosahedron of blobs compared to analytic expressions for the 
Gibbs-Boltzmann distribution. These figures were created using data from 6 independent runs of 4 • 10 5 
time steps with a small time step size of At = 0.04 (6nr]a'j)/ksT to avoid rejections; no rejections occurred 
during these simulations. (Left) Equilibrium distribution of the height h , the distance from the center of 
the icosahedron to the wall. (Right) Equilibrium distribution of the angle 9, where 9 = 0 indicates that the 
heavy blob is at the bottom of the icosahedron, and 9 = n indicates that it is at the top. As expected, we 
see the distribution skewed towards smaller values of 9 due to the gravitational force. 


on the orientation of the sphere. Furthermore, by symmetry, the gravitational force (perpendicular 
to the wall) cannot induce rotation of the sphere, and by symmetry, a torque cannot introduce 
vertical displacements. Because of these special symmetries the parallel MSD is linear for all times 


and therefore (44) gives the long-time quasi two-dimensional diffusion coefficient \ 2 D] this has in 
fact been confirmed experimentally with relatively good accuracy for spheres whose center of mass 
is very close to their geometric center [3]. 


E. Colloidal Boomerang 

The authors of reference [1] perform a detailed experimental study of the quasi-two-dimensional 
translational and rotational diffusion of lithographed symmetric right-angle boomerang colloids 
(see the right panel of Fig. [2]) confined between two closely-spaced microscope slides. Subsequently 
this work was extended to asymmetric (L-shaped) right-angle boomerangs [6j as well as non-right- 
angle boomerangs [32]. Some theoretical analysis is also performed assuming that the overdamped 
dynamics of the particles is strictly two-dimensional. Of course, the actual dynamics of the particles 
is three dimensional, and a complete theoretical or numerical analysis of the diffusive dynamics 
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Figure 7: Mean square displacements for a sphere with nonuniform mass distribution diffusing near a 
planar boundary. Symbols show numerical results obtained by processing long equilibrium trajectories, 
while lines show theoretical predictions. These figures were generated using data from 16 independent 
trajectories of length 5 TO 5 time steps, using the RDF temporal integrator with a small time step size of At = 
0.004 (6nr]a^)/kBT to eliminate rejections. (Left) Parallel (.D||(t)) and perpendicular (D±(t)) translational 
MSD of a rigid icosahedron of blobs. The solid black line shows the theoretical parallel translational MSD 


predicted by (44) for a rigid sphere using the best-known approximations to the hydrodynamic mobility 
flD2|D3|D4| ), while the dashed black line shows the asymptotic perpendicular translational MSD ( [45]) . As 
expected, the icosahedron behaves like a sphere with equal effective hydrodynamic radius [IS]. (Right) 
Parallel (x—x or y—y) as well as perpendicular (z—z) components of the rotational mean square displacement 


(34). The dashed line shows the asymptotic rotational MSD ( |47[ ). We see that the rotational dynamics of 
the rigid icosahedron and a true sphere are also in good agreement. 


requires the complete formalism developed in this paper. 

In this section we examine a single symmetric right-angle boomerang near a single no-slip 
boundary (bottom wall) in the presence of gravity. We choose to study a single boundary rather 
than a slit channel as done in the experiments in order to simplify the hydrodynamic calculations 
of mobilities [35|; in principle one can construct tabulated approximations of self and pairwise 
mobilities in a slit channel but this is quite complex and expensive [62] . While we cannot make 
direct comparisons with the experimental values reported in Ref. [TJ in this work, we can still 
address the fundamental questions about differences between fully three-dimensional and quasi two- 
dimensional diffusion. Specifically, by enlarging the gravitational force we apply to the boomerang 
(i.e., increasing its effective density mismatch with the solvent), we can cause the motion to be 
more or less confined to a two dimensional plane parallel to the bottom wall. In this section we 
use microns as the unit of length, seconds as the unit of time, and milligrams as the unit of mass. 
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For hydrodynamic calculations, we construct a blob model of a boomerang and try to match 
the physical parameters in the experiments [TJ as close as possible. Our model of the boomerang 
particle is constructed by rigidly connecting 15 blobs, one at the cross point, and 7 for each arm, as 
illustrated in the right panel in Fig. [lj Prior investigations in the context of the immersed boundary 
method [56], which we have also confirmed independently by using the Rotne-Prager tensor as the 
pairwise blob mobility, have shown that to construct a good hydrodynamic approximation of a rigid 
cylinder of radius r using blobs, one should set the effective hydrodynamic radius of each blob to 
a ~ 3/2 r, and place the blobs centers on a line at a distance of around a (the precise value does 

not matter much). Following these recommendations, we set the blob radius to a = 0.325, which 
gives an effective cylinder radius of 0.265, and the blobs are spaced a distance 0.3 apart. Note that 
in this minimally-resolved blob model the cross-section of the arms of the boomerang is cylindrical 
rather than square, as would be more realistic for modeling the lithographed particles. We have, 
however, compared to a more resolved 120-blob model constructed from the initial boomerang 
by replacing each of its 15 blobs by 8 smaller blobs of radius 0.1625 placed at the vertices of 
a rectangular prism of size 0.15 x 0.285 x 0.245 centered at the location of the original blob. 
We find only minor differences with the minimally-resolved model, for example, in bulk (without 
confinement) the diffusion coefficients in the plane of the boomerang are computed to be (in units 
of nm 2 /s) 0.243 and 0.283 for the 15-blob model, and 0.245 and 0.291 for the 120-blob model. 

For a free boomerang far away from boundaries, there is a unique CoM that, due to symmetry, 
must lie on the the line that bisects the boomerang arms. Also, there must be a unique point on the 
bisector for which there is no coupling between torque applied out of the plane of the boomerang 
and the translational motion in the plane of the boomerang. We can consider this point as the 
CoH for quasi-two-dimensional diffusion [I], although, as already explained, this point is not a CoH 
in the strict sense for three-dimensional diffusion. The locations of the bulk CoM and the bulk 
quasi-two-dimensional CoH, which we shall henceforth imprecisely refer to as just CoM and CoH, 


can be computed from (40) and (41), respectively. For our blob model, we compute the CoH to 
be is about 1.08 microns away from the cross point (center of the intersection blob), and the CoM 
is 0.96 microns from the cross point; we get the same estimates from the more refined 120-blob 
models. These numbers compare favorably to the experimental findings in [T], where the CoH is 
estimated to be a distance of 1.16 microns from the cross point; the CoM is not mentioned in the 
experimental works on boomerang particles. The difference between the CoH and CoM is too small 
for this specific particle shape for us to be able to tell the difference to within statistical errors; 
in future work we will look for other planar particle shapes for which the difference may be more 
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Figure 8: Translational MSD of a right-angle symmetric boomerang for gravity g = 1 (left panel) and 
g = 20 (right panel) parallel and perpendicular to the wall, see legend. The same trajectories are used but 
the parallel MSD is computed using three different choices of the tracking point: cross point (CP) at the 
corner of the right angle, (one of the) tip(s) of the boomerang, and the center of hydrodynamic stress (CoH), 
which we note is essentially indistinguishable for this purpose from the center of mobility (CoM). For the 
perpendicular MSD we only show the results for the tip and show the expected asymptotic value with a 
dotted line. The slopes of the parallel MSD at the origin are shown with lines (see legend) and estimated 


using Monte Carlo sampling via (44). We see that by calculating the MSD tracking the CoH/CoM, the MSD 


is essentially linear with time, and therefore the short time diffusion coefficient (44) matches the (unique) 
long time diffusion coefficient X 2 D- This is not the case when tracking the cross point or the tip of one of 
the arms. 


significant and measurable in both simulations and experiments. 

The total gravitational force applied to the body is 0.18 x g ( ksT/a ) where g is a parameter 
that we vary; we split the gravitational force evenly among the 15 blobs. Here g = 1 gives a rough 
approximation of the gravitational binding experienced by the actual lithographed particles, which 


have a density of 1.2g/cm 3 . Each blob is also repelled from the wall using the potential (49) with 


screening length b = 0.5a and strength e = 23.08A:T. The gravitational height (46) for one of the 
two (equivalent) tips of the boomerang are shown in Table [T] for several values of g. Since the tips 
are the points that are most likely to venture further from the wall, these values give an indication 
of how close to two-dimensional the dynamics of the boomerang is. 
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9 

X2 D(gm 2 /s) 

X2D/X3D 

Xe (rad 2 /s) 

hg(gm) 

X2D/X3D sphere 

1 

0.226 

0.834 

1.42 

1.77 

0.837 

10 

0.194 

0.716 

0.79 

0.605 

0.724 

20 

0.185 

0.683 

0.22 

0.419 

0.680 


Table I: Long-time quasi two-dimensional diffusion coefficient \ 2 d for the boomerang colloid at g = 1, 10, 


and 20 estimated using (44) and tracking the CoH. The diffusion coefficient for a free boomerang in an 


unbounded fluid \ 3 D is computed using (43). The rotational diffusion coefficient Xo is calculated using 


(Cl). Comparing the effective gravitational heights (|46j), calculated for the tip of the boomerang, and the 
boomerang’s arm length (2.1 gm) or the blob radius (0.325 gm) gives an indication of how flat the boomerang 
is against the wall. In the last column we estimate the reduction in mobility relative to bulk for a sphere of 
the same radius as the blobs and at the same gravitational height as the tip of the boomerang. 


1. Translational Diffusion 

In Fig. [8] we show the parallel and perpendicular translational MSDs of the boomerang for a 
weak (g = 1, left panel) and a strong gravitational sedimentation (g = 20, right panel), where 
strong here means that the gravity is sufficiently strong to keep the boomerang essentially flat 
against the surface. For the parallel MSD, we show results based on three different choices of the 
tracking point: 1) the CoH, or in this example, equivalently the CoM; 2) the center of the blob 
at the tip of one of the arms; and 3) the cross point (CP), which is the center of the blob at the 
cross point where the arms meet. For an unconfined boomerang, we expect that the parallel MSD 
measured using the CoM will be strictly linear in time. We also expect that for the boomerang 
confined to a plane by strong gravity and exhibiting quasi two-dimensional diffusion, the MSD 
will be linear in time when tracking the CoH mi- However, in our simulations, we find that 
due to the close proximity of the CoH and the CoM, the MSD is identical to within statistical 
error independent of which of these two points was tracked. For clarity, we only include the MSD 
calculated from the CoH in our results, with the understanding that the CoM is indistinguishable 
at this level of accuracy. We see from the figure that by choosing the CoH as the tracking point, we 
obtain a MSD that is linear over all times up to statistical accuracy for both gravities. This means 
that we can get an accurate estimate of the long-time diffusion coefficient X 2 D by using equation 


(44) over a broad range of gravities. This statement should, however, be checked for other particle 


shapes for which the CoH and CoM are sufficiently far apart, before drawing broad conclusions. 
In Table |T] we show the estimated long-time parallel diffusion coefficient X 2 D obtained from 


(44) for different strengths of the gravitational sedimentation. We find that, perhaps surprisingly, 
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Figure 9: Planar rotational mean squared displacement (501 of the colloidal boomerang for three different 
strengths of the gravitational confinement. 


the presence of the boundary does not strongly reduce the effective short-time diffusion coefficient 
compared to bulk, except at the largest gravity. In the last column of the table we give the 
corresponding reduction in mobility for a sphere of the same radius as the blob radius, as obtained 


from the theoretical estimate (D4) averaged against a Gibbs-Boltzmann distribution ~ exp (—h/h g ) 
for the height above the wall; a remarkable agreement is observed despite the significant difference 
in the particle shape. The value of the quasi two dimensional diffusion coefficient is measured 
experimentally in Ref. [1] for the case of a boomerang particle confined between two microscope 
slides a distance 2 fim apart, and a value of X 21 F = 0.054 fim 2 /s is reported. This is lower 
than the values calculated here, which we expect is largely due to the absence of the top wall 
in our simulations, which will significantly increase the drag on the boomerang for such strong 
confinement. 


2. Rotational Diffusion 


To estimate the quasi-two-dimensional rotational diffusion coefficient xe-, measured experimen¬ 
tally in Ref. [1], we project the bisector of the boomerang arms into the x — y plane, and define 
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9(t) to be the angle of rotation around the z axis between this projected vector at time t and the 
initial projected bisector. We count each full counterclockwise rotation as an addition of 2ir, and 
similarly we subtract 2ir for each full clockwise rotation, allowing the value of 9 to take values in 
all of M. Truly two-dimensional rotational diffusion (where the colloid stays in the x — y plane) 
corresponds to 9(t ) being standard Brownian motion. We define a planar rotational mean square 
displacement from increments of the angle 6, 


D e (' T) = (M 2 ) = [(9(t + T)-9(t)f). 


(50) 


In Fig. [9] we show numerical results for Dg (r). We see that for g = 20, when the boomerang is 
most flat, the quantity Dg is linear in time for all times to within statistical error bars, while for 
lower gravities we see some deviations from linearity, as expected since the definition of 9 assumes 
the diffusion is essentially two-dimensional. 

We estimated the short-time planar rotation coefficient 


Xe 


lim 

At->0 


(A 9 2 ) 
2A t 


using Monte Carlo averaging based on (Cl) from Appendix |C| and tabulate the computed values in 
Table [IJ Note, however, that the trajectory Dg (r) is only continuous if the boomerang never flips, 
i.e., the bisector is never nearly perpendicular to the wall; note that we do observe flips for two lower 
gravities. We see that xe is much larger for lower gravities, both because the boomerang diffuses 
more rapidly far from the wall, and also because in low gravity, the boomerang is not confined to the 
x — y plane, and hence small changes in orientation can lead to large changes in our calculated two- 
dimensional angular displacement. The rotational diffusion coefficient measured experimentally in 
Ref. [1] for a boomerang confined between two microscope glass slips is xe = 0.044rad 2 /s, which 
is much smaller than our result for g = 20. This is most likely in large part due to the absence of 
drag from the top wall. Additionally, without this second boundary, our simulated boomerang is 
able to rotate out of the x — y plane more easily, reaching configurations where a small change in 
orientation can lead to a large change in 9(t). 


V. CONCLUSION 

In this paper, we studied the Brownian motion of rigid bodies of arbitrary shape immersed in 
a viscous fluid in the overdamped regime, in the presence of confinement and gravity. We pa¬ 
rameterized the orientation of the rigid bodies with normalized quaternions, which offer several 
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advantages over other previously-used representations. Furthermore, we do not assume any par¬ 
ticular symmetry for the rigid bodies, and we account carefully for the fact that the hydrodynamic 
mobility N depends on the configuration due to confinement or hydrodynamic interactions with 
other particles. We derived the appropriate form of the overdamped Langevin equations of motion, 
including all of the stochastic drift terms required to give the correct Gibbs-Boltzmann distribution 
in equilibrium, and to preserve the unit norm constraint of the quaternions. 


In section III we developed temporal integrators for the rigid-body overdamped Langevin system 
and presented two ways to handle the stochastic drift term. The first approach is a generalization of 
the well known midpoint Fixman scheme [sami, which generates the drift terms using a midpoint 
predictor step but requires a costly application or factorization of N -1 . The second approach 
employs a Random Finite Difference approach to generate the drift terms using only applications 
of N and N 2 , making it an appealing choice. The RFD approach is especially promising in 
situations where the action of the mobility and the stochastic terms are generated by using a 
fluctuating hydrodynamics fluid solver, as in the fluctuating force coupling method (FCM) [4T], or 
extensions of our fluctuating immersed boundary (FIB) method jl5j| to include rotlet (and possible 
also stresslet) terms in the minimally-resolved blob model. 

In Section |TV| we performed several numerical simulations of the Brownian motion of rigid par¬ 
ticles diffusing near a wall in the presence of gravity, motivated by a number of recent experiments 
studying the diffusion of asymmetric spheres [3], clusters of spheres M, and boomerang colloids 
mi- First, we examined the behavior of a tetramer formed by rigidly connecting four colloidal 
spheres together, modeling colloidal clusters that have been manufactured in the lab ESj. Sec¬ 
ond, we studied the rotational and translational diffusion of a colloidal sphere with nonuniform 
density, modeling recently-manufactured “colloidal surfers” [8j in which a dense hematite cube is 
embedded in a polymeric spherical particle. Finally, we investigated the quasi two-dimensional 
diffusive motion of a dense boomerang colloid sedimented near a no-slip boundary, inspired by 
recent experimental studies of lithographed boomerang-shaped particles mm- 

We demonstrated that the choice of tracking point is crucial when computing the translational 
diffusion coefficient, as already observed and explained in Refs. HI El- In particular, we demonstrate 
that in some cases there exists a suitable choice of the origin (around which torques are expressed) 
which can be used to obtain an approximate but relatively accurate formula for the effective long¬ 
time diffusion coefficient in the directions parallel to the boundary. For highly symmetric shapes 
with a clear geometric center it turned out that the “obvious” tracking point is the best one to 
use. However, for the boomerang shapes studied here we found that the CoH and CoM are so 
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close to each other that we cannot numerically distinguish between them. Therefore, it remains to 
be confirmed whether the CoH, rather than the CoM, is the correct point to track in quasi-two- 
dimensional confinement as claimed in Refs. Bi- Ideally one would find a planar shape for which 
these two points are far apart; unfortunately our calculations indicate that all of the boomerang 
shapes studied in published experiments have a CoH and a CoM that are too close to each other 
to be distinguished to within experimental and statistical accuracy. Additional investigations of 
other particle shapes are necessary to reach more definitive conclusions about diffusion in quasi- 
two-dimensional (strong) confinement. 

In many practical situations only part of diffusing particle may be tracked, for example, a unit 
of a protein may be labeled by a fluorescent dye. In such cases, one must be very careful in 
interpreting the results for translational diffusion as if the particle were spherical and the center 
of the sphere were tracked. Furthermore, there are many particle shapes for which there is no 
obvious geometric center and it is then not trivial to determine what the best point to track is, 
even if one can track an arbitrary point on the body. In general, we find that there is no exact 
closed-form expression for the long-time quasi-two-dimensional coefficient; it appears necessary to 
perform numerical simulations in order to study the long-time diffusive dynamics of even a single 
rigid body in the presence of confinement. Our temporal integrators can easily be extended to study 
quasi-two-dimensional suspensions of passive or active particles sedimented near a boundary, which 
is quite relevant in practice since active particles often have metallic components and are therefore 
much denser than the solvent M. 

In our simulations, the time step size was strongly restricted in order to keep the rigid body 
from passing through the wall. To this end, we rejected steps that encountered an unphysical 
state (e.g., a configuration where the computed mobility matrix is not positive semi-definite). This 
naive approach modifies the dynamics in a way that violates ergodicity and detailed balance, and 
we reduced our time step size to avoid performing a significant number of rejections. Several 
more sophisticated approaches exist that may solve this problem, including Metropolization m , 
adaptive time-stepping |T6], or continious-time discretizations [63]. Employing these techniques in 
our integrators remains an area of future exploration. 

Recently, the Brownian motion of a spheroid (an axisymmetric particle) near a single no-slip 
wall has been studied El by using a finite element method for pre-computing the hydrodynamic 
mobility over many positions of the particle relative to the wall, and using the RFD approach to 
compute the divergence of the mobility in expectation. The strategy of Ref. El of pre-computing 
the mobility does not extend to suspensions of particles, and constructing body-fitted finite element 


43 


meshes and solving the resulting Stokes equations is rather computationally intensive. In this work 
we relied on a simple rigid multiblob approach for computing the hydrodynamic mobilities (52] . 
using direct dense linear algebra to compute inverses and Cholesky factorizations. This was useful 
for validating our methods, but it does not scale well with increasing numbers of rigid bodies or blobs 
per rigid body. Furthermore, the analytical approximation we used for the blob mobility is valid 
only for the case of a single no slip boundary [55], and even in that case it is not guaranteed to lead 
to a symmetric positive-definite grand mobility for all configurations. The RDF scheme developed 
in this work can be coupled with a computational fluid solver, similarly to the approach taken in 
the FIB method in a way that will allow us to do simulations in more complex geometries 
such as slit or square channels or chambers, and scale to large numbers of blobs. The required 
rigid-body immersed boundary method has recently been developed m , and in the future the 
temporal integrators developed in this work will be employed to account for the Brownian motion 
of the rigid particles. 
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Appendix 

Appendix A: Quaternions and Rotation 

In this appendix we derive some relations regarding the quaternion representation of orienta¬ 
tions, as used in the main text. 


1. Rotating a body 


In this section, we derive eq. (10). We proceed by first writing the Rotate procedure ([9]) using 
its definition, and then expand the trigonometric functions to second order. Letting 6 = {s,p}, 
and ||cl7|| = co, we have 


Rotate(0, uAt) = 


s cos (^^) - p ■ sin <jj/w 
s sin (^j^) uj/uo + cos (^j^) p + sin (^^) uj X p/u 


s( 1- ^ 
+ (l - 


2 At 2 ' 


p-uf 


2 At 2 


At 


p + UJ X P±f 


s - p ■ u>^ 


UJ- 


'■At 2 


p + suf - 

_ T . (u> ■ cj) At 2 

=6 + '&uAt - - --- 


s 

P 


+ O (At 3 ) 
+ O (At 3 ) 


e + o (At 3 ). 


2. Torques 

In this section, we consider the case when a torque is generated by a conservative potential 
U v (p), so that r = —dU^/dp. Here p represents the oriented angle associated with orientation. 
For the purposes of this discussion, we neglect the dependence of potential on location, since this 
will have no bearing on the torque. Consider extending the energy to depend on a quaternion U (6) 
such that when ||0|| = 1, we have U v (p) = U{6^). We want to be able to write the torque, r in 
terms of U(6) without needing to convert first to p. 

Only quaternions with unit norm represent a viable orientation, and therefore the value of the 
potential off of this constraint has no physical meaning and should not affect the torque in any 
way. The projected gradient of U (6) on the unit 4-sphere is 

B K = d K_(f) d E\p = P d R 

dG dG V 96) Pe dG ’ 
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where Pg = I — 66 1 , and it can easily be checked that HU Pg = 'I'' . In section II A, we saw that 


T 

~P 

dp = 7) 

—p ■ dp 

Sl + P 

Z 

sdp +p x dp 


dO = - 
2 


so that the change in potential energy due to a small rotation dip is 


du = -T-d v = 9 dL-de = -\ 


dll , dU . , , , 

lte P - d( P- -^-(sdV + Pxdv) 


Using the vector identity a ■ (b X c) = c • (a x 6) , we can rewrite this as 


T.<V=- 


leading to the identification of torque as 


dU , r dU 
— P -(S -P)^ 


■ dp, 


1 

T ~ 2 


dU , r 91 / 

—p-(s - p )q£ 


MU 


dU 


MU 


= -ty 1 — = -Hr 

50 0 06 06 


(Al) 


Appendix B: Stochastic Drift Terms 


Here we show that the temporal integrators introduced in Section III generate the correct 
stochastic drift terms, more precisely, they are first-order weakly accurate integrators. We will find 
it convenient in the following calculations to consider the drift term separated into multiple pieces 


as done in (28). We first derive (13), which we use in the following subsections. We start with the 


form of the drift written in (12), denoting for simplicity M = and using indicial notation 


with Einstein’s implied summation convention for clarity, 


0 e ■ M 


= Oj yMijj = 0 j (^ ik M k i^ji) 

= (OjVik) + 4 >ik ( djM kl ) + Vv k M kl (OjVji) 

= (dj'&ik) Mki'bji + ^ik (djMkl) 'J 'ji 
= - ~gM kk 6i + ( OjM k i ) 4 >ji 

where we used ([8]) to go from the second to the third line. To go from the third to the fourth 
line we used the relationship (Oj'&ik) 4 >ji = —5 k idi/A, which can be shown by a straightforward 
calculation. In (somewhat ambiguous) matrix notation, we can write 


0 e ■ M = T' (OqM) : - -Tr ( M ) 0, 


(Bl) 


which we use in proving first-order weak accuracy of our numerical schemes next. 
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We also derive a similar relation for the drift including translational degrees of freedom, as 


given in (28). We use Einstein’s implicit summation notation, where Greek indices range over 
components of location q , s,t , and u range over components of orientation 6 , and i represents any 
component of x. We now expand the z-tli component of the stochastic drift using the chain rule 
and (l8b , 




f)+d s {M^ st ) 
d a + d s (* it M“JV su 

(dpM-f) + {d s MY) *st 

'bis (d a M?f) + (dsMSn*, 


+ 


— ( dnN m p ) Z. n p A 




(d s V it )(M“T)* s 


(B2) 


1. Fixman’s Method 


To show that (26) is equivalent to (12), we can use the general identity that given two matrices 
A ( x ) and B (x), 


A o B W (d x A) : ( BB t A t ) + AB W 

= ^ (d x • ( ABB t A t ) - A d x ■ ( BB t A t )) + AB W, 


(B3) 


in law, where {(S^A) : (BB 1 A 1 )} . = ( diAij) BjkB m kAi m , and we used the product rule to obtain 


the second line of (B3). Applying this identity we obtain 


y/2k B TENoN-2W 

= (k B T) 8 X ■ (SNS t ) - (k B T) HAT {d x -S T ) + y/2k B TEN*W 
= ( k B T) 8 X • (SNS t ) + y/2k B TEN*W, 


(B4) 


where we used 


To show that scheme (27) produces the correct drift terms, we consider the drift for q and 6 


separately. In the following expression, Greek indices range only over components corresponding 
to q and not those corresponding to 6 , and indices s, t, u, and v correspond to only components 
of 6 . All other indices range over every variable. Letting Ax P,n+ ^ = y/k B TAt ZkiNlWZ 1 be 
the stochastic term from the increment x P ^ n+ ^ — x^, the stochastic drift generated for q by the 
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corrector stage is equal to 

Ath<£ =Vk B TAt(d k N aj ) A a% n+ *N~* (w?' 1 + IF/ 1 ’ 2 ) 

=fc s TAt [id,N a] ) Nj m W.Z 1 + (9 s iV ai ) VstN^WgA x 

where all matrices are evaluated at x n and the term involving dUi comes from expanding the 
Rotate procedure in the predictor stage. After taking expectation and noting that, for example, 
N a p = , we obtain the stochastic drift 

(A t hC) = ksTAt {{dpMZf) + (d a M%) **) , (B5) 


which matches the first row of the second line in (B2) as required. 


For the drift in the 6 direction, we expand the Rotate procedure in the corrector stage to first 
order in At to obtain 


n, 2 
3 


A th 9»"=*.i ( yk„TU ( a t N „) Axf ” + =JVy ’ (W"’ l + W. 

_ keTAt + w „,2^ N i N i ^ w n.l + w n,2 


7 - 71,1 


= (k B TAt) * st \^d a N ti ) N^Wj; 

+ ( d u N ti ) VuvN^WZ’ 1 ) N~^ (iF”’ 1 + W n ' 2 


k B T At 


3 1 3 

wJ 1 ’ 1 + iyjivj (if ;- 1 + w”’ 2 


e,. 


After taking expectation we obtain the deterministic drift 
<A th 0”) = k B TAt 


** (d a M?f) + vk st * u „ - 


(B6) 


(B7) 


which matches the second row of the second line in (B2) as required. Note that a direct application 
of the Euler-Heun scheme [52] to the (26) would require the final update of orientation to be 


e n +1 = e n + ^P’ n+ lu>P’ n+ lAt. 


Using \I>P ’ n+ 2 instead of \E ,ri here generates the drift term —Tr (iVf) 9/ 4; here we obtain that part 
of the stochastic drift by using the Rotate procedure instead of a simple additive update of the 
quaternions. 
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2. Random Finite Difference Scheme 

To show that the random finite difference term generates the correct drift, we need to show that 
<5 -1 ^N — N n ^j W is a good approximation to d x ( N ) : S T in expectation. We use the convention 
that Greek indices correspond to translational degrees of freedom, s and t correspond to angular 
degrees of freedom, and the remaining indices are summed over all variables. Expanding the RFD 
term gives 

A th x? = ^ fa - Wj =^j-d k (Nij) Ax k Wj + 0(6), 
where Ax k = x k — x^. Expanding the increment Ax k gives 

(Nij) Ax k Wj =k B T (d a (Nij) W a Wj + d s (N t] ) T st W t W ’■) . 

Taking expectation gives the desired result 

<A th s?> =k B Td k (Nij)Z kj . 


Appendix C: Planar Rotational Diffusion Coefficient 


In Section IV E we computed the two dimensional rotational diffusion coefficient by measuring 
the change in angle 6 of the bisector of the boomerang projected onto the x — y plane. This is 
a convenient notion of rotational diffusion when the boomerang lies flat, in which case it can be 
used to measure the z — z component of However, in the case of general three dimensional 

motion, the relationship between Dg(r) and is not as simple. In this appendix, we derive 

the relationship between the short time diffusion coefficient in 6 and the rotational mobility. 

We consider the boomerang at an initial location and orientation, x = (q,s,p) and let v = 
(vi,V 2 ,vs) be the unit vector pointing in the direction of the bisector. We let M = M LOT (x) be 
the rotational mobility evaluated at the initial configuration. Finally, we let Q be the projection 
operator that projects vectors onto the xy plane. We then consider the change in the angle 8 
between the projected bisector and the x axis after a rotation over a small time increment At. Let 
<f> be the angle of rotation over this time increment, and R be the rotation matrix that applies this 
small rotation. For small 0, the change in the scalar angle 8 is 

|| (Qv) x (QRv) 


A 8 =1 


\Qv\ 


+ o 
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For small <f>, we approximate the rotation matrix as 


R= 2 


pp T + sP + [ s 2 - ^ ) I 


=i+$+o 


where $ is the cross product matrix for 0, i.e.,<h;E = <f> x x. Using this approximation to R, we 
get an expression for the instantaneous planar diffusion coefficient xe ( x ), 

(A.0 2 ) 

Jim . 2 k B TAt =Ms3 _ “ (^ 3 Mi 3 + 2v 2 v 3 M 23 ) + a 2 v 2 (vjM u + v\M 22 + 2Mi 2 ViU 2 ) , (Cl) 

where a = (v 2 + v 2 ) 1 . The average short time projected rotational diffusion coefficient is then 
X 0 = (xo ( x)), where the average is taken with x distributed according to the equilibrium Gibbs- 


Boltzmann distribution (24) 


Appendix D: Hydrodynamic Mobility of a Sphere Near a Wall 


A low-order approximation of the perpendicular and parallel translational mobilities of a sphere 
next to a no-slip boundary is derived by Swan and Brady (45] as a generalization of the Rotne-Prager 
tensor using Blake’s image construction [64] . Neglecting stresslet contributions, this approximation 
gives the self-mobility 


p±(h) _ 9a a 3 a 


Po 

M|| (h) 


1 8 h + 2/i 3 8 h 5 

9 a 2a 3 o 
= 1 - — + 


(Dl) 


(jlq 16 h 16 h 3 16/i 5 ’ 

where po = (67rr/a) _1 is the mobility in an unbounded domain. We do not reproduce the lengthier 
formula for the other components of the mobility. 

More accurate formulas for the self-mobility of a sphere near a wall are available. A very 
good approximation to the perpendicular mobility is given by a semi-empirical rational relation 
approximation to an exact series of Brenner (65] , 

M-lW 6 (f) + 2 (a) / D2) 

^ ^ (a) 2 + 9 (a) + 2 ' 

The hard sphere approximation to the parallel mobility is given by a combination of a near-wall 
expression derived using lubrication theory and a truncated expansion in powers of a/h which is 
more accurate further from the wall. The near-wall calculation given by Goldman and Brenner 
gives 

M|| (h) 2 (In (£) -0.9543) 


Mo 


(In (S)) 2 - S.lSSln (|) + 1.591 


(D3) 
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and it is used when h — a < 0.03a. When the sphere is further from the wall, we calculate the 
parallel mobility from the exact power series expansion truncated to fifth order mi 

Ho 16 h 8 h 3 256 h 4 16 h 5 

We were unable to find more accurate expansions for the rotation-rotation and rotation- 
translation components of the mobility of a sphere near a wall. Therefore, we compute them 
from a cubic spline fit to the numerical mobility obtained from a sphere discretized with 162 blobs. 
We observe that this well-resolved multiblob model provides a rather accurate approximation, as 
confirmed by comparing the numerical translational mobilities to the above theoretical expansions. 
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Note that, in principle, the formalism developed here can directly be applied to two dimensions by 
replacing quaternions with complex numbers; a rotation of <fi radians in a counterclockwise direction is 
associated with the complex number 6$ = exp (i(p) = cos <t> + i sin cf>. 

If the accumulation of numerical errors has caused|||0" +1 1| — l| > e, for some tolerance e, one should 
renormalize the quaternion, 6 n+1 6 n+1 / ||0 71+1 ||. 

Numerically, a uniformly-distributed unit 4-vector can be sampled by generating a vector of 4 standard 
Gaussian random variables and normalizing the result; to see this observe that the resulting distribution 
must be uniform by virtue of the rotational invariance of the multivariate Gaussian distribution. 

The Euler-Heun method is the natural generalization of the Euler-Maruyuama method to SDEs with 
Stratonovich noise [12]. 

For three dimensional particles diffusing parallel to the x — y plane, one can define a quasi-two- 
dimensional CoH as an origin for which a torque around the z axis does not induce any translational 
motion in the x — y plane; the set of such points is a line parameterized by a parameter s, 


r CoH (s) = (-( M Uz f - M UtTx s)/M UzT „, (M Ux p x + Mu z r s)/M ul ^, s ). 


The two-dimensional result (411 is a special case of this more general formula for s = 0. 



